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ABSTRACT 

We present galaxy stellar mass functions (GSMFs) at 2 = 4-8 from a rest-frame ultraviolet (UV) 
selected sample of ~4500 galaxies, found via photometric redshifts over an area of ~280 arcmin 2 in 
the CANDELS/GOODS fields and the Hubble Ultra Deep Field. The deepest Spitzer/ IRAC data 
yet-to-date and the relatively large volume allow us to place a better constraint at both the low- 
and high-mass ends of the GSMFs compared to previous space-based studies from pre-CANDELS 
observations. Supplemented by a stacking analysis, we find a linear correlation between the rest- 
frame UV absolute magnitude at 1500 A (Muv) and logarithmic stellar mass (logAf*) that holds for 
galaxies with log (M^/Mq) < 10. We use simulations to validate our method of measuring the slope 
of the log Af*-Afuv relation, finding that the bias is minimized with a hybrid technique combining 
photometry of individual bright galaxies with stacked photometry for faint galaxies. The resultant 
measured slopes do not significantly evolve over 2 = 4-8, while the normalization of the trend exhibits 
a weak evolution toward lower masses at higher redshift. We combine the logAf*-Afuv distribution 
with observed rest-frame UV luminosity functions at each redshift to derive the GSMFs, finding that 
the low-mass-end slope becomes steeper with increasing redshift from a = —1.551q q® at 2 = 4 to 
a = —2.25^0 35 at 2 = 8. The inferred stellar mass density, when integrated over Af* = 10 8 -10 13 
Mq, increases by a factor of 10^2° between 2 = 7 and 2 = 4 and is in good agreement with the time 
integral of the cosmic star formation rate density. 

Subject headings: galaxies: evolution — galaxies: formation — galaxies: high-redshift — galaxies: 
mass function 


1. INTRODUCTION 

The near-infrared (near-IR) capability of the Wide 
Field Camera 3 (WFC3) on board the Hubble Space 
Telescope ( HST ) has now yielded a statistically signifi¬ 
cant sample of galaxies in the early universe, enabling 
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us to pass the era of simply discovering very distant 
galaxies and enter an era where we can perform system¬ 
atic studies to probe the underlying physical processes. 
Such studies have begun to make progress toward under¬ 
standing galaxy evolution at high redshift, with particu¬ 
lar advances in measurements of the rest-frame ultravio¬ 
let (UV) luminosity function (e.g., Bouwens et al. 2011; 
Oesch et al. 2012; Schenker et al. 2013b; Lorenzoni et 
al. 2013; Finkelstein et al. 2015; Bouwens et al. 2015), 
the rest-frame UV spectral slope (e.g., Finkelstein et al. 
2012; Bouwens et al. 2014), and the cosmic star forma¬ 
tion rate density (SFRD; see Madau & Dickinson 2014 
for a review). 

A key constraint on galaxy evolution that has only re¬ 
cently begun to be robustly explored is that of the growth 
of stellar mass in the universe. This measurement re¬ 
quires a combination of deep rest-frame UV data with 
constraints at rest-frame optical wavelengths, to better 
probe the emission from older, lower-mass stars. The 
mass assembly history across cosmic time is governed by 
complicated processes, including star formation, merg¬ 
ing of galaxies, supernova feedback, etc. In spite of the 
complexity of the baryonic physics of galaxy formation, 
however, various studies have found that global galaxy 
properties, such as star formation rate (SFR), metallic- 
ity, size, etc., all correlate tightly with the stellar mass 
(Kauffmann et al. 2003b; Noeske et al. 2007; Tremonti et 
al. 2004; Williams et al. 2010), implying that the stellar 
mass plays a major role in galaxy evolution. Thus, de- 
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termining the comoving number density of galaxies in a 
wide range of stellar masses (i.e., the galaxy stellar mass 
function [GSMF]) and following the evolution with red- 
shift constitutes a basic and crucial constraint on galaxy 
formation models. Specifically, obtaining robust obser¬ 
vational constraints on the low-mass-end slope of the 
GSMF can provide insights on the impact of feedback 
on stellar mass buildup of low-mass galaxies (Lu et al. 
2014), as current theoretical models predict steeper low- 
mass-end slopes than those previously observed (e.g., Vo- 
gelsberger et al. 2013; see Somerville & Dave 2014 for a 
review). Consequently, the evolution of GSMFs over the 
past 11 billion years has been extensively investigated ob- 
servationally during the past decade (e.g., Marchesini et 
al. 2009; Baldry et al. 2012; Ilbert et al. 2013; Muzzin et 
al. 2013; Tomczak et al. 2014), mostly using traditional 
techniques (e.g., 1/V m ax, maximal likelihood) that assess 
and correct for the incompleteness in mass. 

Nonetheless, the GSMF in the first 2 billion years after 
the Big Bang has remained poorly constrained, due to (i) 
limited sample sizes, particularly of low-mass galaxies at 
high redshift; (ii) systematic uncertainties in stellar mass 
estimations; and (iii) the lack of Spitzer Space Telescope 
Infrared Array Camera (IRAC; Fazio et al. 2004) data 
with comparable depth to HST/ WFC3. IRAC data are 
essential to probe stellar masses of galaxies at z > 4 as 
the 4000 A/Balmer break moves into the mid-infrared 
(mid-IR), beyond the reach of HST (and ground-based 
telescopes). 

An alternative approach to deriving the GSMF takes 
advantage of the fact that the selection effects and in¬ 
completeness are relatively well known and corrected for 
in rest-frame UV luminosity functions. Therefore, one 
can derive GSMFs by taking a UV luminosity function 
and convolving it with a stellar mass versus UV lumi¬ 
nosity distribution at each redshift. Using this approach, 
Gonzalez et al. (2011) derived GSMFs at 4 < z < 7 of 
Lyman break galaxies (LBGs) over ~33 arcmin 2 of the 
Early Release Science (ERS) field (for z = 4-6, and the 
Hubble Ultra Deep Field [HUDF] and the Great Obser¬ 
vatories Origins Deep Survey [GOODS] fields with pre- 
WFC3 data for z = 7) utilizing the HST/ WFC3 and 
IRAC GOODS-South data. They reported a shallow 
low-mass-end slope of —(1.4-1.6) at z = 4-7, but the 
small sample size and limited dynamic range made it 
difficult for them to explore the mass-to-light distribu¬ 
tion at z > 4, and their GSMFs were derived under an 
assumption that the mass-to-light distribution at z ~ 4 
is valid up to z 

More recently, several studies have utilized the large 
HST data set from the Cosmic Assembly Near-infrared 
Deep Extragalactic Legacy Survey (CANDELS; Grogin 
et al. 2011; Koekemoer et al. 2011) to make progress 
on the measurements of the GSMFs at z > 4. Dun¬ 
can et al. (2014) and Grazian et al. (2015) derived the 
GSMFs of galaxies at 4 < z < 7 in the GOODS-S 
field and GOODS-S/UDS fields, respectively, using the 
CANDELS HST data and the Spitzer Extended Deep 
Survey (SEDS; Ashby et al. 2013) IRAC data. These 
studies have reported a steeper low-mass-end slope of 
a ~ —(1.6-2.0) than previous studies, but the uncer¬ 
tainties are still large, and the inferred evolution of the 
low-mass-end slope of the GSMFs remains uncertain. 


Here we probe galaxy buildup from z = 4 out to 
z = 8, aiming to improve on the limiting factors dis¬ 
cussed above, with a goal of providing robust constraints 
on the GSMFs of galaxies at 4 < z < 8. We do this 
by combining near-IR data from CANDELS GOODS- 
S and GOODS-N fields with the deepest existing IRAC 
data over the GOODS fields from the Spiteer-CANDELS 
(S-CANDELS; PI Fazio; Ashby et al. 2015) and the 
IRAC Ultra Deep Field 2010 Survey (UDF10; Labbe et 
al. 2013). We obtain reliable photometry on these deep 
IRAC data by performing point-spread function (PSF) 
matched deblending photometry, enabling us to extend 
the exploration of the GSMFs to lower stellar masses and 
higher redshifts than previous studies. Furthermore, a 
special emphasis is put on quantifying and minimizing 
the systematics inherent in our analysis via mock galaxy 
simulations. While taking a similar approach of convolv¬ 
ing a rest-frame UV luminosity function with stellar mass 
to rest-frame UV luminosity distribution as Gonzalez et 
al. (2011), the increased sample size and deeper data en¬ 
able us to bypass the limitations of the previous studies, 
as we measure the mass-to-light ratio distribution at ev¬ 
ery redshift. 

This paper is organized as follows. Section 2 intro¬ 
duces the HST data sets used in this study, as well as 
our sample at 3.5 < z < 8.5 selected by photometric 
redshifts, and discusses our IRAC photometry, which is 
critical for the stellar mass estimation described in Sec¬ 
tion 3. Section 4 presents stellar mass versus observed 
rest-frame absolute UV magnitude (M*-Muv) distribu¬ 
tions and introduces a stacking analysis and mock galaxy 
simulations. By combining the rest-frame UV luminos¬ 
ity function with the M*-Afuv distribution, we derive 
GSMFs and stellar mass densities in Section 5 and 6, 
respectively. A discussion and summary of our results 
follow in Section 7 and Section 8. Throughout the pa¬ 
per, we use the AB magnitude system (Oke & Gunn 
1983) and a Salpeter (1955) initial mass function (IMF) 
between 0.1 Af© and 100 M©. All quoted uncertainties 
represent 68% confidence intervals unless otherwise spec¬ 
ified. We adopt a concordance ACDM cosmology with 
Hq = 70 = 100/r km s -1 Mpc -1 , Om = 0.3, and 12 a = 
0.7. The HST bands F435W, F606W, F775W, F814W, 
F850LP, F098M, F105W, F125W, F140W, and F160W 
will be referred to as B, V, i, Igu, z, 1098; Y> J> JHuq, 
and H, respectively. 

2. DATA 

Constraining GSMFs requires a deep multiwavelength 
data set over a wide area in order to probe the full dy¬ 
namic range of a galaxy population. In this section, we 
describe the HST imaging used to select our galaxy can¬ 
didates, as well as the candidate selection process. We 
then discuss the procedures used to measure accurate 
photometry for these galaxies from the Spitzer /IRAC S- 
CANDELS imaging. 

2.1. HST Data and Sample Selection 

The galaxy sample employed in this study is from 
Finkelstein et al. (2015), to which we refer the reader 
for full details of the HST data used and the galaxy sam¬ 
ple selection. This sample consists of ~7000 galaxies 
selected via photometric redshifts over a redshift range 
of 2 ; = 3.5-8.5. These galaxies were selected using HST 
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Fig. 1. — Example of our IRAC photometry modeling procedure. From left to right. (1) //-baud WFC3 imaging of a l'x 1' region in the 
GOODS-S field, (2) S-CANDELS 3.6 fi m imaging, (3) the best-fit 3.6 fi m model image, and (4) the T-PHOT residual image (i.e., real science 
image subtracted by the best-fit model image). Our high-redshift galaxies (gray squares) are often blended with nearby foreground sources; 
therefore, we perform PSF-matched photometry on the S-CANDELS data using the WFC3 H -band imaging as a prior on the position 
and morphology of sources. Using the T-PHOT software package, we convolve the if-band image with empirically derived IRAC PSFs to 
generate low-resolution (IRAC) model images, allowing the flux of each source to vary to simultaneously fit all sources in the IRAC data. 


imaging data from the CANDELS (Grogin et al. 2011; 
Koekemoer et al. 2011) over the GOODS (Giavalisco 
et al. 2004) North (GOODS-N) and South (GOODS-S) 
fields, the ERS (Windhorst et al. 2011) field, and the 
HUDF (Beckwith et al. 2006; Bouwens et al. 2010; Ellis 
et al. 2013) and its two parallel fields (Oesch et al. 2007; 
Bouwens et al. 2011). 14 We use the full data set which in¬ 
corporates all earlier imaging from HST Advanced Cam¬ 
era for Surveys (ACS), including the B, V, i, Isi 4 , and 
z filters. We also use imaging from the HST WFC3 in 
the Yogs, F. J, Tffi 4 o, and H filters. A complete de¬ 
scription of these data is presented by Koekemoer et 
al. (2011, 2013). These combined data are three-layered 
in depth, comprising the extremely deep HUDF, moder¬ 
ately deep CANDELS-Deep fields, and relatively shallow 
CANDELS-Wide and ERS fields, designed to efficiently 
sample both bright and faint galaxies at high redshifts. 

As described by Finkelstein et al. (2015), sources were 
detected in a summed J + H image, using a more aggre- 
sive detection scheme compared to that used to build the 
official CANDELS catalog (e.g., Guo et al. 2013, Barro 
et al. 2016, in preparation) to detect faint sources at 
high redshifts, yielding a catalog with the total number 
of sources a factor of two larger. 15 To minimize the pres¬ 
ence of spurious sources, only those objects with > 3.5cr 
significance in both the J and H bands were evaluated 
as possible high-redshift galaxies. Photometric redshifts 
were estimated with EAZY (Brammer et al. 2008), and 
selection criteria devised based on the full redshift prob¬ 
ability density function (pdf; P(z)) were applied (Finkel¬ 
stein et al. 2015). These criteria were designed to assess 
the robustnes of the sample, requiring the primary red- 
shift solution to have more than 70% of the integrated 
redshift pdf, and the narrowness of P{z ), by limiting the 
sample to those for which the integral of the redshift pdf 
for the corresponding redshift bin is at least 25% of the 
total integral of the pdf. A comparison of our photomet- 

14 Finkelstein et al. (2015) also include about 500 more galaxies 
from the Abell 2744 and MACS J0416.1-2403 parallel fields from 
the Hubble Frontier Fields data set, while we do not, reducing our 
full sample to —7000 galaxies. 

15 Specifically, Finkelstein et al. (2015) used DE- 
TECT.MINAREA = 7 and DETECT.THRESH = 0.6, while 
the “hot” mode of the official CANDELS catalog was built with 
DETECT_MINAREA = 10 and DETECT.THRESH = 0.7. 


ric redshifts for 171 galaxies with spectroscopic redshifts 
in our 2 > 3.5 sample shows an excellent agreement, with 
cr{tSz/(l + ^spec)) = 0.03 (after 3cr clipping). All of the 
selected sources were visually inspected for removal of 
artifacts and stellar contaminants. Active galactic nuclei 
(AGNs) identified in X-rays were also excluded from the 
sample. Our final parent sample consists of 4156, 2056, 
669, 284, and 77 galaxies at z = 4 (3.5 < z < 4.5), 5 
(4.5 < z < 5.5), 6 (5.5 < z < 6.5), 7 (6.5 < z < 7.5), and 
8 (7.5 < z < 8.5), respectively. 

2.2. IRAC Data and Photometry 

At 3.5 < z < 8.5, the observed mid-IR probes rest- 
frame optical wavelengths redward of the Balmer/4000 A 
break. Deep Spitzer/ IRAC data are therefore critical to 
constrain stellar masses and the resulting GSMF. One 
of the key advances of our study is the significantly in¬ 
creased depth in the 3.6 and 4.5 pm IRAC bands pro¬ 
vided by the new S-CANDELS survey (Ashby et al. 
2015). The final S-CANDELS mosaics in the GOODS- 
S and GOODS-N fields (where the former includes the 
HUDF parallel fields) include data from three previous 
studies: GOODS, with integration time of 23-46 hr per 
pointing (Dickinson et al., in preparation); a 5 ; x 5' region 
in the ERS observed to 100 hr depth (PI Fazio); and the 
IRAC Ultra Deep Field 2010 program, which observed 
the HUDF and its two parallel fields to 120, 50, and 100 
hr, respectively (Labbe et al. 2013). The total integra¬ 
tion time within the S-CANDELS fields more than dou¬ 
bles the integration time for most of the area used in this 
study (to >50 hr total), reaching a total formal depth of 
26.5 mag (3<r) at 3.6 pm and 4.5 pm (Ashby et al. 2015). 
Imaging at 5.8 and 8.0 pm was obtained with IRAC as 
part of the GOODS program. However, these data are 
too shallow to provide meaningful constraints for high- 
redshift faint galaxies, so we do not include them in our 
analysis. 

The full-width at half-maximum (FWHM) of the PSF 
of the IRAC data (~1 , . , 7 at 3.6 pm versus ~0'/19 at 1.6 
pm with WFC3) results in non-negligable source confu¬ 
sion, making accurate flux determinations challenging, 
as shown in Figure 1. The second panel of Figure 1 
shows that our high-redshift galaxies are extremely faint 
and are often blended with nearby bright sources in the 
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mid-IR. making simple aperture photometry unreliable. 
For reliable IRAC photometry on the deep S-CANDELS 
data, we therefore perform PSF-matched photometry us¬ 
ing the T-PHOT software (Merlin et al. 2015), an updated 
version of TFIT (Laidler et al. 2007), on the S-CANDELS 
3.6 and 4.5 /xm mosaics. This PSF-matched photom¬ 
etry uses information in a high-resolution image (here 
the H band), such as position and morphology, as priors. 
Specifically, we use isophotes and light profiles from the 
detection ( J+H) image obtained by the Source Extractor 
package (SExtractor; Bertin & Arnouts 1996). The high- 
resolution image was convolved with a transfer kernel to 
generate model images for the low-resolution data (here 
the IRAC imaging), allowing the flux in each source to 
vary. This model image was in turn fitted to the real 
low-resolution image. The IRAC fluxes of sources are 
determined by the model that best represents the real 
data. 

As the PSF FWHM of the high-resolution image (H 
band) is negligible (~(y/19) when compared to those of 
the low-resolution IRAC images (~l , /7), we used IRAC 
PSFs as transfer kernels. We derive empirical PSFs in 
each band and in each field by stacking isolated and mod¬ 
erately bright ([3.6], [4.5] = 15.5-19.0 mag) stars identi¬ 
fied in a half-light radius versus magnitude diagram in 
IRAC imaging. As t-phot requires the kernel to be in 
the pixel scale of the high-resolution image, each star 
image was 10 times oversampled to generate the final 
PSFs in the same pixel scale of the 77-band image (0 / /06 
pixel -1 ). They were then registered, normalized, and 
median-combined to generate the final IRAC PSFs. 

Preparatory work for running T-PHOT includes per¬ 
forming background subtraction on the S-CANDELS 
mosaics using the script subtract_bkgd.py (written by 
H. Ferguson; see Galametz et al. 2013 for details) and 
dilation of the J+H SExtractor segmentation map, anal¬ 
ogous to the traditional aperture correction. The need 
for this dilation step originates from the fact that un¬ 
der nonzero background fluctuations, isophotes of faint 
sources include a smaller fraction of their total flux com¬ 
pared to bright sources, and therefore their IRAC fluxes 
based on isophotes tend to be underestimated. To in¬ 
clude the light in the wings and to counterbalance the 
underestimation of flux for faint sources of which the 
amount depends on the isophotal area, an empirical cor¬ 
rection factor to enlarge the SExtractor segmentation 
map was devised by the CANDELS team by compar¬ 
ing isophotal area from deep and shallow images (see 
Galametz et al. 2013 and Guo et al. 2013 for details). 
We applied this empirical correction factor to the J+H 
segmentaion map using the program dilate (De Santis 
et al. 2007) while preventing merging between sources. 

To correct for potential small spatial distortions or 
mis-registrations between the high-resolution and low- 
resolution images, a second run of t-phot was performed 
using a shifted kernel built by cross-correlating the model 
and real low-resolution images. Figure 1 presents an ex¬ 
ample of our IRAC photometry procedure on a subregion 
in GOODS-S. With the exception of very bright sources, 
the residual image is remarkably clean, highlighting the 
accuracy of this procedure. 

2.2.1. Verification 
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Fig. 2. — Comparison between our T-PHOT S-CANDELS pho¬ 
tometry and the official CANDELS TFIT SEDS photometry from 
Guo et al. (2013) and Barro et al. (2016, in preparation) for IRAC 
3.6 /mi (upper) and 4.5 /mi (lower) bands for sources with S/N 
> 1 (black), S/N > 3 (green), and S/N > 5 (red) in both cata¬ 
logs. Blue circles and error bars indicate the median and robust 
standard deviation of the magnitude difference in each magnitude 
bin for sources with S/N > 1, and the blue dashed lines encom¬ 
pass the central 68% of the distribution. The red vertical dotted 
line denotes the 5er S-CANDELS depth. The bias seen at faint 
magnitudes is due to Eddington bias of upscattered sources in the 
shallower SEDS data, which is shown as the gray shaded area (for 
the S/N > 1 cut; see Section 2.2.1 for more details). The SEDS 
data are shallower; thus, the agreement between the two catalogs 
in magnitude bins brighter than this range indicates that our pho¬ 
tometry is accurate. 


We tested the accuracy of our IR AC photometry in two 
ways. First, we compared our catalog with the official 
CANDELS catalogs (Guo et al. 2013; Barro et al. 2016, 
in preparation) in which IRAC fluxes were obtained from 
the shallower SEDS data using TFIT. Figure 2 shows the 
comparison between our magnitude and that from the of¬ 
ficial catalog for sources with signal-to-noise ratio (S/N) 
greater than 1, 3, and 5 in both catalogs. As we compare 
two catalogs generated from images with different depths 
with a certain S/N cut, faint sources are dominated by 
background fluctuations and the Eddington bias of up- 
scattered sources in the shallower SEDS data, making 
the comparison unreliable (see Figure 13 of Guo et al. 
for a similar trend and discussion of the flux compari¬ 
son). Taking a similar approach to that of Guo et al. 
(2013), we estimated this magnitude range in which the 
flux comparison is unreliable due to the Eddington bias 
for the S/N > 1 cut. We first found the best-fit power law 
to the differential number count density of the sources 
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in the official CANDELS catalog without any S/N cut. 
Then, we compared the differential number count den¬ 
sity of sources with S/N > 1 with the best-fit power law 
to find the magnitude that the former starts falling be¬ 
low 80% of what is expected by the best-fit power law 
([3.6] = 25.1 and [4.5] = 25.5; gray shaded region in 
Figure 2). In magnitude bins brighter than this range, 
the comparison indicates an excellent agreement with a 
negligible systematic offset. 

Second, we performed a mock source simulation in or¬ 
der to validate our photometry. Briefly, mock sources 
with varying physical properties (e.g., size, light profile, 
luminosity, redshift) were generated using the GALFIT 
software (Peng et al. 2002), of which flux densities in each 
band are assigned using the updated Bruzual and Char¬ 
iot (CB07; Bruzual & Chariot 2003) stellar population 
synthesis (SPS) models. While the exact shape of the as¬ 
sumed distribution of physical parameters such as size or 
light profile can impact the detection rate of sources close 
to the sensitivity limit and thus the results of simulations 
designed to correct for completeness, t-phot photome¬ 
try is by design limited to the sources recovered in the 
high-resolution detection image. Therefore, our simula¬ 
tion results should not be sensitive to those assumptions 
in the first order. These mock sources were convolved 
with the PSF of each filter and added at random loca¬ 
tions in the (77-band PSF-convolved) high-resolution and 
IRAC low-resolution images. Our simulation thus ac¬ 
counts fully for source confusion with nearby foreground 
real sources, but we constrain the number density of 
mock sources to be negligible (5 arcmin -2 ) compared to 
that of real sources to ensure that our simulation results 
are not dominated by self-crowding among the inserted 
artificial sources. The mock sources were then recovered 
using the same procedures as for real sources, including 
T-PHOT photometry. Figure 3 shows a comparison be¬ 
tween the input and recovered 3.6 /mi and 4.5 /i.m mag¬ 
nitudes. It is encouraging that we find no systematic 
offset between the two down to the 50% completeness 
limit of 25 mag (Ashby et al. 2015), given that the IRAC 
photometry is affected with many sources of uncertainty 
that demand accurate background subtraction, aperture 
correction scheme (dilation), and buildup of transfer ker¬ 
nels. As any observed offsets from these simulations are 
not statistically significant, we do not apply any correc¬ 
tion to our observed IRAC photometry. 

2.2.2. Visual inspection 

A practical limitation exists when building empirical 
IRAC PSFs for deep fields, as such surveys by design 
target a relatively small area well off the Galactic plane, 
resulting in few stars present in the imaging. The se¬ 
vere source confusion in the IRAC data further reduces 
the number of isolated stars that can be used for the 
creation of PSFs. Therefore, although our PSF-nratched 
IRAC photometry significantly improves photometric ac¬ 
curacy over more traditional methods (Lee et al. 2012), 
our IRAC photometry may be imperfect. Its significance, 
however, is likely small, as already discussed in the pre¬ 
vious sections. 

Another source of uncertainty is that TFIT/t-phot as¬ 
sumes no morphological ^-corrections (no variation in the 
surface profile or morphology) between short- and long- 
wavelength images, which is likely not the case (although 



Fig. 3.— Comparison between the input and recovered IRAC 
3.6 /im {upper) and 4.5 pm {lower) magnitude in our mock source 
simulations. Symbols are color-coded by blendedness in the two 
filters (shown in the inset in binary notation)—00: contamination- 
free; 01: contaminated in 4.5 /am; 10: contaminated in 3.6 /im; 11: 
contaminated in both 3.6 and 4.5 pm. Large black circles and error 
bars represent the median and standard deviation of the magnitude 
difference of sources in each magnitude bin with the magnitude 
difference between the input and the recovered less than 2 mag, 
demonstrating the reliability of our IRAC photometry down to the 
50% completeness limit of 25 mag (red vertical dotted line; Ashby 
et al. 2015). 

we attempted to mitigate this by using the reddest HST 
image as our high-resolution image). However, our high- 
redshift galaxies are small enough to not be resolved in 
the low-resolution image (see Figure 25 of Ashby et al. 
2013), and thus it should have a negligible effect on the 
derived fluxes. Bright (and extended) foreground sources 
in close proximity to our high-redshift sample, however, 
are prone to imperfect modeling in this case (even with 
perfect PSFs) and leave residuals that in turn can signif¬ 
icantly affect the photometry of faint sources nearby. 

To account for these uncertainties, we visually in¬ 
spected the IRAC science and residual images of all 
~7000 sources in our sample to ensure that their IRAC 
photometry is reliable. Sources falling on strong residu¬ 
als from nearby bright sources often have recovered S/Ns 
that are too high (even when we cannot visually identify 
counterparts in the IRAC images) or significantly nega¬ 
tive, which indicates that their photometry is not reliable 
in general. This is confirmed by a mock source simula¬ 
tion in which we added mock sources at various posi¬ 
tions around a bright source, finding that the recovered 
flux is highly biased (either overestimated or underesti¬ 
mated depending on the “yin” and “yang” of the residual 
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on which the mock source was inserted). We therefore 
caution against blindly taking IRAC fluxes from a cat¬ 
alog and stress the importance of visual inspection of 
IRAC images and residuals of all objects, as contami¬ 
nated sources can significantly impact studies on indi¬ 
vidual galaxies or with small number statistics (e.g., the 
high-mass end of the GSMF). 

Contaminated sources are excluded from the sample 
in our subsequent analysis. This leaves ~63%, 63%, 
54%, 61%, and 57% of our z = 4,5, 6, 7,8 parent sample 
free from a possible contamination from nearby bright 
sources, 16 resulting in 2611, 1292, 364, 172, and 44 
sources in our final z = 4, 5,6, 7,8 selection. Among 
the final sample, 1172/2611, 480/1292, 108/364, 41/172, 
and 6/44 (45%, 37%, 30%, 24%, and 14%) sources at 
z = 4,5,6, 7,8, respectively, have > 2cr detections at 3.6 
//m, and 613/2611, 211/1292, 43/364, 12/172, and 1/44 
(23%, 16%, 12%, 7%, and 2%) show S/N > 5. A final 
multiwavelength catalog was constructed by combining 
the HST catalog and the IRAC t-phot catalog for this 
final sample. 

3. STELLAR POPULATION MODELING 

We derived stellar masses and rest-frame UV luminosi¬ 
ties for our sample galaxies by fitting the observed spec¬ 
tral energy distribution (SED) from the B, V, i, I% 5 o, 
z, y 0 98 , Y, J, J7/i4o, H, 3.6 /im, and /.5 /im data to 
the Bruzual & Chariot (2003) SPS models. We refer the 
reader to Finkelstein et al. (2012, 2015) for a detailed ex¬ 
planation of our modeling process. Briefly, we modeled 
star formation histories (SFHs) as exponentially declin¬ 
ing (t = 1 Myr-10 Gyr), constant (r = 100 Gyr), and 
rising (r = -300 Myr, -1 Gyr, -10 Gyr). Allowable 
ages spanned a lower age limit of 1 Myr to the age of 
the universe at the redshift of a galaxy, spaced serni- 
logarithmically, and metallicity ranged from 0.02 to 1 
Z 0 . We assumed the Calzetti et al. (2000) attenuation 
curve with E(B — V) = 0.0-0.8 (Ay = 0.0-3.2), and 
intergalactic medium (IGM) attenuation was applied us¬ 
ing the Madau (1995) prescription. All the models were 
normalized to a total mass of 1 Mq . 

One of the major sources of uncertainty in stellar 
mass measurements of high-redshift galaxies derived 
from broadband imaging data via SED fitting is the con¬ 
tribution of nebular emission. Spectroscopically measur¬ 
ing the contribution of strong emission lines (e.g., Ha, 
[O ill]AA4959, 5007) to the broadband fluxes for z > 4 
galaxies is not currently feasible because they redshift 
into the mid-IR. Many efforts, however, have been made 
taking an alternative approach of investigating IR AC col¬ 
ors of spectroscopically confirmed galaxies at the redshift 
range such that only one of the first two IR AC bands is 
expected to be contaminated by a strong emission line 
(e.g., Shim et al. 2011; Stark et al. 2013). Together with 
an inference from the direct measurements of nebular 
lines of galaxies at lower redshift (e.g., z ~ 3; Schenker 

16 These uncontaminated source fractions are consistent with 
the confusion-free fraction of the S-CANDELS images found by 
Ashby et al. (2015). While the symbols in Figure 3 are color-coded 
by the blendedness, which is determined based on the SExtractor 
segmentation map, in this section we visually inspected individ¬ 
ual sources to exclude sources for which photometry is affected by 
residuals of nearby bright sources. That is, those excluded by our 
visual inspection are a subset of sources that are color-coded as 
being contaminated in Figure 3, and likely catastrophic outliers. 
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Fig. 4.— Comparison of the inferred rest-frame 

EW([0 ill]AA4959, 5007) distrbution derived from our SED 

fitting analysis for galaxies in our 2 : ~ 4 sample with spectroscopic 
redshifts (black solid histogram) with the spectroscopic EW([0 ill]) 
distribution of 20 LBGs from Schenker et al. (2013a, blue dashed 
histogram; eight galaxies out of their 28 targets were not de¬ 
tected). The good agreement implies that our implementation of 
nebular emission lines in our stellar population models is a fair 
representation of reality. 


- This study 

.Schenker+13 (3.0 < z < 3.8) 



et al. 2013a), several studies claim that the contribution 
of nebular emission to the broadband fluxes for galaxies 
at high redshift, and thus to the inferred stellar mass, 
can be significant (e.g., Schaerer & de Barros 2010). 

While a more concrete answer on the nebular contribu¬ 
tion will need to wait for the advent of the James Webb 
Space Telescope (JWST ), we took into account the con¬ 
tribution of nebular emission in our SED modeling in a 
self-consistent way following the prescription of Salmon 
et al. (2015). In this prescription, the strengths of H 
and He recombination lines (including H/l) are set by 
the number of non-escaping ionizing photons, and the 
strengths of metal lines relative to H/3 are given by In- 
oue (2011). The nebular line strengths are thus a func¬ 
tion of the population age and metallicity (which sets 
the number of ionizing photons), as well as the ioniz¬ 
ing photon escape fraction. We added the nebular lines 
to the stellar continuum assuming an escape fraction of 
zero 17 and no extra dust attenuation for nebular emis¬ 
sion (i.e., E(B - V) stellar = E(B - F) nebular)- Fig¬ 
ure 4 presents a comparison of the inferred equivalent 
width (EW) distribution of [O ill]AA4959, 5007 from our 
SED modeling for spectroscopically confirmed galaxies at 
3.5 < z < 4.5 with the observed EW([0 ill]) histogram of 
20 LBGs at similar redshifts (3.0 < z < 3.8) in Schenker 
et al. (2013a), showing a good agreement. As Figure 4 is 
based on spectroscopically confirmed galaxies, they are 
biased toward high-mass (median ± standard deviation 
of the logarithmic stellar mass = 9.7±0.5) and UV-bright 
(—21.3 ± 0.8) galaxies. However, the two samples shown 
in Figure 4 have similar stellar mass distributions. 

Our final stellar+nebular line stellar population mod- 

17 Many studies at z ~ 3, the highest redshift where we can 
estimate the LyC escape fraction before the IGM becomes opaque 
to ionizing photons, have shown no evidence of high escape fraction, 
only placing an upper limit of <0.1 (lcr; Siana et al. 2015). 



















The Galaxy Stellar Mass Function at z = 4-8 


7 


els are integrated through all of the filter bandpasses in 
our photometry catalog. The best-fit model for each 
source was found as the one that best represents the 
observed photometry via y 2 minimization. During this 
procedure, we accounted for uncertainties in the zero- 
point and aperture corrections by adding 5% of the flux 
in quadrature to the flux error in each band. For our 
fiducial stellar mass for each galaxy, we adopted the me¬ 
dian mass obtained from a Bayesian likelihood analysis 
following Kauffmann et al. (2003a). We describe the pro¬ 
cedure briefly here, but we refer the reader to Kauffmann 
et al. (2003a) and our previous work (Song et al. 2014) 
for more details of the analysis (see also Salmon et al. 
2015; Tanaka 2015). In short, we used the y 2 array that 
samples the full model parameter space of our SPS mod¬ 
els to compute the four-dimensional posterior pdf of free 
parameters (dust extinction, age, metallicy, and SFH) 
using the likelihood of each model, £ oc e~ x / 2 . We as¬ 
sumed flat priors in parameter grids and z = z p hot- The 
stellar mass for each grid point is taken to be the normal¬ 
ization factor between the observed SED and the model. 
Then, the 1-dimensional posterior pdf for stellar mass 
was obtained by marginalizing over all the parameters. 
The median and the central 68% confidence interval of 
stellar mass were computed from this marginalized pdf. 

The rest-frame absolute magnitude at 1500 A, M\jy, 
was obtained from the mean continuum flux density of 
the best-fit model in a AA res t = 100 A band centered 
at rest-frame 1500 A. Its uncertainty was derived from 
100 Monte Carlo simulations in which the redshift un¬ 
certainty was accounted for by varing the redshift in our 
Monte Carlo simulations following the pdf, P(z), ob¬ 
tained from our photometric redshift analysis (Finkel- 
stein et al. 2015). Systematic biases and uncertainties of 
our SED fitting method were estimated via mock galaxy 
simulations and will be discussed in Section 4.3. 

4. STELLAR MASS-REST-FRAME UV LUMINOSITY 
DISTRIBUTION 

With the stellar mass (A/*) and the rest-frame absolute 
UV magnitude (Afuv) of our galaxy sample measured 
from the previous section, we now explore the correlation 
between these properties in our sample to infer whether 
it is possible to derive a scaling relation. 

4.1. Overall M* -Afuv Distribution 

Figure 5 presents the stellar mass versus rest-frame UV 
absolute magnitude distribution at each redshift. Over¬ 
all, we find a strong trend between stellar mass and rest- 
frame absolute UV luminosity at all redshifts probed in 
this study. The scatter (standard deviation) in logarith¬ 
mic stellar mass is about 0.4 dex (0.52, 0.42, 0.36, 0.40, 
and 0.30 dex at z = 4, 5, 6, 7, and 8, respectively, mea¬ 
sured as the mean of standard deviation in logarithmic 
stellar mass in bins with more than five galaxies), and no 
noticeable correlation of the scatter is found with red¬ 
shift or UV luminosity. The scatter at the bright end 
(measured at —21.5 < Afuv < —20.5), where the effect 
of observational uncertainty should be minimal, is 0.43, 
0.47, 0.36, 0.52, and 0.40 dex at z = 4, 5, 6, 7, and 8, 
respectively, similar to the quoted value above and the 
scatter at the faint end (at —19.0 < Afuv < —18.0) of 
0.51, 0.39, 0.39, and 0.41 dex at z = 4, 5, 6, and 7, 
respectively. 


Often raised as a weakness in studies of rest-frame 
UV-selected galaxies is that such studies by construc¬ 
tion miss dusty star-forming or quiescent galaxies. As 
our sample is selected mainly from rest-frame UV, the 
observed A/*-A/uv distribution shown in Figure 5 is also 
susceptible to this weakness. Interestingly, however, we 
do observe populations of UV-faint galaxies with high 
mass (given their UV luminosity) on the upper right 
part of the Af*-A-fuv plane at z = 4 and 5. Although 
the lower limit of SFRs inferred from the UV luminos¬ 
ity and the Kennicutt (1998) conversion assuming no 
dust (upper x-axis in Figure 5) indicates that they are 
not completely quenched systems, their inferred (dust- 
uncorrected) SFRs are down to 2 orders of magnitude 
lower than those on the A-f^-Afuv relation (to be derived 
in the following section). Outliers off the best-fit relation 
by more than 1 dex make up 6% (155/2624) of the total 
sample at z = 4. The fraction decreases as redshift in¬ 
creases to 3% (12/365) at z = 6. This increasing fraction 
of massive and UV-faint galaxy populations from z = 6 
to z = 4 implies either that we may be witnessing the 
formation of dusty star-forming or quiescent populations 
that are very rare at high redshift (z ~ 6-7) or that the 
duty cycle of those populations at high redshift is lower 
than that at low redshift (with a star-forming time scale 
much longer than 100 Myr) so that fewer such galaxies 
are observed with the current flux limit at higher redshift. 
Given the young age of the universe at high redshift, the 
latter requires a very early and fast growth of stellar 
mass for those galaxies that are completely quenched or 
highly dust extincted by z ~ 7. As we do not see such 
extremely UV-luminous populations at higher redshifts 
in the UV luminosity functions (e.g., Finkelstein et al. 
2015; Bouwens et al. 2015), we regard the former as a 
more plausible scenario. 

Interestingly, we do not find such populations in the 
opposite (low-mass) side at a given UV luminosity. The 
lack of bright and low-mass galaxies in the lower left re¬ 
gion of Figure 5 is most clearly shown at z = 4. This is 
unlikely to be a selection effect or observational uncer¬ 
tainty, as had there been galaxies with log(Af* /Mq) > 9, 
they should have been detected in both WFC3/IR and 
IRAC; although we impose an S/N cut in both J and 
H bands in our sample selection and may thus be bi¬ 
ased against the bluest galaxies, this only applies for UV 
bins fainter than those discussed here. It should not be 
an artifact of our SPS modeling, as the minimum mass- 
to-light ratio allowed in our SPS models is well below 
the mass-to-light ratio distribution of our sample. Lee 
et al. (2012), based on LBGs selected at z = 4-5 over 
the GOODS field, interpreted the absence of undermas- 
sive galaxies as evidence of smooth growth for UV-bright 
galaxies that has lasted at least a few hundred million 
years. If dust extinction is proportional to UV luminos¬ 
ity (e.g., Bouwens et al. 2014), this indicates a lack of 
UV-bright galaxies with very high specific SFRs (sSFR 
= SFR/Af*). This lack of UV-bright and low-mass galax¬ 
ies at all redshifts we probe provides a further support 
on our claim in the paragraph above of a growing pop¬ 
ulation of dusty star-forming or quiescent galaxies seen 
between z = 6 and z = 4. 

4.2. Stacking Analysis 
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Fig. 5.— From upper left to lower right , stellar mass versus rest-frame UV absolute magnitude at 1500 A (Mjjv) at z = 4—8. For 
reference, SFR inferred from the UV luminosity and the Kennicutt (1998) conversion assuming no dust is shown in the upper x-axis. Small 
gray filled circles indicate objects with IRAC detection (> 2cr at 3.6 pm), while gray open circles are those with nondetections in IRAC 
(< 2cr at 3.6 pm). Gray error bars represent the 68% confidence intervals in stellar mass and rest-frame absolute UV magnitude. Large 
red filled circles are the median stellar masses in each rest-frame absolute UV magnitude bin of 0.5 mag, while large open circles indicate 
bins containing a single galaxy. Black error bars are standard deviations in stellar mass in each UV luminosity bin. Blue stars indicate 
median stellar masses in each rest-frame absolute UV magnitude bin from our median-flux stacking analysis in Section 4.2, with error bars 
denoting the lcr uncertainty, including both photometric error and sample variance. We derived the best-fit relation (red solid line) by 
fitting data points that combine red filled circles with blue filled stars in a redshift-dependent UV magnitude range specified in Section 4.3 
(indicated as the light-red and light-blue filled regions). The lcr uncertainty of the best-fit M*— Myy relation is denoted as the light-red 
shaded region. The gray arrows and horizontal error bars at the bottom show the characteristic UV magnitude, L* , of the UV luminosity 
function (Finkelstein et al. 2015) at each redshift. Gray dotted lines indicate minimum mass-to-light ratio allowed in our SPS models. We 
find that the best-fit relation has a nonevolving slope at z = 4—6, which is marginally steeper than a constant mass-to-light ratio (gray 
dashed line; normalized to the mass-to-light ratio of the Milky Way), and shows a weak evolution in the normalization. The inferred stellar 
mass from the best-fit M*-Muy relation for galaxies with M\jy = —21 (~ L*) is log(M*/M 0 ) = 9.70, 9.59, 9.53, 9.36, and 9.00 at 2 : = 4, 
5, 6, 7, and 8, respectively. The typical mass-to-light ratio of galaxies at 2 : = 4—8 at the rest-frame UV absolute magnitude of the Milky 
Way (Muv = —20.5) is lower by a factor of ~30 (130) at 2 : = 4 (8) than that of the Milky Way. 
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Fig. 6 .— Median flux-stacked SEDs at z = 4-8 in bins of rest-frame absolute UV magnitude with AMuy = 0.5 mag, for bins with 
M\jy < —17 and more than 10 (for z = 4—5) or 5 (for z = 6-8) galaxies (corresponding to blue stars in Figure 5). The stacked SEDs are 
denoted by filled circles and downward-pointing arrows (indicating 2cr upper limits for bands with S/N < 2), with the rest-frame absolute 
UV magnitudes given by the inset text. The solid lines and open squares indicate the best-fit SPS models and model bandpass-averaged 
fluxes, respectively. The best-fit SPS models for stacked SEDs with S/N < 2 in 3.6 ^m are shown as dotted lines, indicating that the 
inferred mass for stacked SEDs with S/N < 2 in 3.6 /im is uncertain. 


Despite our deep IRAC data, individual galaxies in 
our sample, especially those in faint UV luminosity bins 
(which likely have low masses), often suffer from low S/N 
in the IR AC mosaics. This can make the reliability of an 
M*-Muv relation derived based on individual galaxies 
questionable. We therefore performed a stacking analy¬ 
sis to increase the S/N and to examine the typical stellar 
mass in each rest-frame absolute UV magnitude bin. We 
built median flux-stacked SEDs, comprising a total of 12 
bands (B, V, i, 7 8 i4, 2 , Vigs, Y, J, JH 140 , H, 3.6 /im, 4.5 
/jm), for galaxies in each UV magnitude bin of 0.5 mag in 
the full sample. Uncertainties on the stacked SEDs were 
assigned as the quadrature sum of the photometric uncer¬ 
tainty and the uncertainty due to sample variance (het¬ 
erogeneity of the SEDs of galaxies) estimated via boot¬ 
strapping on galaxies in each UV magnitude bin. The 
latter dominates the error budget at z = 4, contributing 


on average 80% to the total uncertainty, but the contri¬ 
bution of sample variance decreases with redshift, down 
to 45% at z = 8. This is a combined effect of (i) de¬ 
creasing outliers in the M*-Muv plane (i.e., decreasing 
fraction of dusty star-forming or quiescent galaxy popu¬ 
lation) with redshift, as seen in Figure 5, and (ii) increas¬ 
ing photometric uncertainty with redshift, as the galaxies 
are fainter compared to the photometric depth. 

Stacked SEDs were analyzed through our SED fitting 
procedures described in Section 3 with the redshift of 
model templates fixed to the median photometric red¬ 
shift of galaxies in each stack. Bands not common to 
all galaxies (e.g., JH 140 covering only the HUDF) were 
included in the SED fitting only when more than half of 
the stacked galaxies have measurements. Our choice is a 
trade-off between making the most of the available infor¬ 
mation and minimizing the chance of biasing our results 
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by including a subset not having the characteristics of the 
parent sample. As the median flux-stacked SEDs short- 
ward of the Lya line may still have some fluxes depending 
on the redshift distribution of the stacked galaxies and 
thus may not represent an SED of a galaxy at the median 
redshift, we excluded bands shortward of the Lya line. 

The best-fit SPS models and stacked SEDs in each red¬ 
shift bin are shown in Figure 6. Overall, the shapes of 
the SEDs are nearly similar but show a weak trend with 
UV luminosity such that the typical UV-bright galaxies 
have slightly redder rest UV-to-optical (observed NIR- 
to-MIR) color than UV-faint galaxies at a given redshift, 
being in qualitative agreement with other previous stud¬ 
ies (e.g., Gonzalez et al. 2012). This trend indicates that 
on average UV-faint galaxies have (mildly) lower mass- 
to-light ratios than UV-bright galaxies, suggesting that 
an Ai*-Afuv relation with a constant mass-to-light ratio 
would not provide a good description of our data. 

The results from our SED fitting analysis on the me¬ 
dian flux-stacked SEDs are included in the _A/*-Muv 
plots of Figure 5. Comparing the stacked points to 
the median of individual galaxies shows that they are 
largely consistent with each other at the bright end 
(Afuv 1$ —(20-21)). But for fainter UV luminosities, 
the stacked points are generally lower than the medians. 
This may reflect the fact that the stellar masses of indi¬ 
vidual IRAC-undetected galaxies are on average biased 
toward higher masses, while they are relatively robust 
for IRAC-detected galaxies. 

4.3. Mock Simulation with Semi-Analytic Models 

As noted in the previous section, stellar mass estima¬ 
tion involves various sources of uncertainty, which im¬ 
pact the derived GSMFs. In our methodology of con¬ 
volving the Af*-Afuv distribution with the completeness- 
corrected UV luminosity functions to derive GSMFs, the 
reliability of our to-be-derived GSMF and its low-mass- 
end slope is tied to our ability to recover the intrinsic 
slope, normalization, and scatter of the Af*-Afuv distri¬ 
bution. 

To explore the systematics and uncertainties in our 
observed M*-Muv distribution introduced by the pho¬ 
tometric uncertainty of our data and our stellar mass 
estimation procedure, and to examine whether we can 
recover the intrinsic Af*-Afuv distribution (and the sub¬ 
sequent GSMF), we simulated Af*-A4uv planes using 
mock galaxies drawn from semianalytic models (SAMs). 
We used synthetic galaxy photometry from the SAMs of 
Somerville et al. (2016, in preparation). These SAMs are 
implemented on halos extracted from the Bolshoi dark 
matter simulation (Klypin et al. 2011), which has a dark 
matter mass resolution of 1.35 xlO 8 h~ 1 Mo and a force 
resolution of 1 A -1 kpc. Galaxies hosted in a halo more 
massive than 10 10 Mq are included in the mock catalog. 
This corresponds to a typical stellar mass of a few times 
10' Mq at z = 4-8 (e.g., Behroozi et al. 2013), similar to 
the minimum stellar mass in our real sample. The most 
notable characteristic of these SAMs is that their light 
cones are specifically designed to provide realizations of 
the hve CANDELS fields (with albeit a factor of 6-9 
larger areas than the actual CANDELS HST coverage), 
aiming to help with interpretation of observational data. 
Moreover, using synthetic galaxy photometry from the 
SAMs has an advantage over using SPS models in that 


mock galaxies have more realistic SEDs based on more 
complicated SFHs and metal enrichment histories, thus 
representing real galaxy populations more closely. We 
refer the reader to Somerville et al. (2012) and Lu et al. 
(2014) for full details of their mock galaxy models. Here 
we specifically use the mock light cone of the CANDELS 
GOODS-S field for our simulation. 

We generated mock galaxy samples by populating the 
Af*-Aiuv plane at each redshift with objects from the 
SAM catalog similar to our real sample in both sam¬ 
ple size and rest-frame UV absolute magnitude distribu¬ 
tion, but with various input slopes ranging from —0.3 
to —0.8. We assumed a lognormal distribution around 
a linear log Af*-Afuv relation with a dispersion of ~0.3 
dex, motivated by the functional form of the observed 
star-forming main sequence at lower redshifts (Speagle 
et al. 2014 and references therein). Although the SAMs 
have an inherent M/L relation, this does not affect our 
results as we randomly draw galaxies from the SAM to 
fill in our simulated plane based on the M/L slope in 
a given simulation. We then added noise in each band 
based on the flux uncertainty of our real sample at a 
given magnitude and perturbed the simulated photome¬ 
try assuming a Gaussian error distribution. The stellar 
masses and rest-frame UV absolute magnitudes of these 
mock galaxies were calculated in the same manner as in 
the real data. The above precedures describe a single 
mock realization of one intrinsic Af*-Afuv distribution. 
For each realization of a given intrinsic A4*-Afuv distri¬ 
bution, we measured the recovered slope, normalization, 
and scatter of the Af*-A/uv distribution. We repeated 
the above procedures 50 times for each input slope and 
redshift, from which we constructed pdfs of the recovered 
slopes and normalizations. 

In addition to allowing us to explore the best way to 
minimize the bias and uncertainty in recovering the in¬ 
trinsic Af*-Afuv distribution, these simulations also al¬ 
low us to test the validity of our stellar mass measure¬ 
ments. First, we find that if we use the classical maxi¬ 
mum likelihood estimator (i.e., the best-fit model) as a 
hducial stellar mass, the uncertainty in the inferred stel¬ 
lar mass for individual galaxies is considerable, differing 
up to 1 dex for galaxies with \og{M* / Mq) ~ 9 at 2 = 4. 
As a result, the scatter of the recovered A/*-Muv dis¬ 
tribution increases significantly compared to the original 
one (a 0.2 dex increase at Afuv ~ “20 and larger for 
fainter galaxies; see also Salmon et al. 2015). This large 
spread in stellar mass in the recovered distribution makes 
it hard to recover the intrinsic relation at 2 > 6 where 
the photometric uncertainty is large and the sample size 
is small. The median mass from our Bayesian likelihood 
analysis described in Section 3 performs much better in 
the sense that it does not significantly increase the scat¬ 
ter of the recovered Af*-Muv plane even in faint UV 
luminosity bins; the scatter of the recovered AL*-A-fuv 
distribution compared to that of the input distribution 
shows an increase of only 0.05-0.10 dex at Afuv ~ —20, 
and the scatter remains nearly constant in fainter UV lu¬ 
minosity bins. However, there exists a noticeable bias in 
the recovered stellar mass for galaxies fainter (and lower 
in mass) than a certain UV magnitude threshold at each 
redshift (hereafter referred to as Afthreshi(-z))- This is be¬ 
cause, for low-S/N data, the stellar mass is determined 
by the assumed priors (e.g., flat priors in parameter grids 
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Fig. 7.— Probability density function of the recovered Af*— 
Muv slope as a function of the intrinsic slope at 2 = 4—8 in our 
SAM mock galaxy simulations. Darker gray colors indicate a higher 
probability. For reference, a 1:1 line is shown as the black solid line. 
The best-fit slope of the Af*—Afuv relation and its 1 a uncertainty 
obtained from our real sample (using our optimized fitting method) 
at each redshift are denoted as the red dashed horizontal line and 
shaded region, respectively. At z = 4—6, we can recover the intrinsic 
slope within ±0.07. At 2 : = 7—8, the uncertainties become much 
larger; thus, as discussed in Section 4.4, we fix the slope (which 
does not significantly evolve from 2 : = 4 to 6 ) at 2 : = 7 and 8 to the 
2 : = 6 value. 

in the case of our SPS modeling) and the data have little 
constraining power. 

Although the recovered M*-Muv distribution above 
Mthreshi (z) is reliable in both bias and scatter, the small 
dynamic range above this limit at high redshift still 
makes the uncertainty of the recovered M* -Muv rela¬ 
tion large. We thus utilize a stacking analysis (described 
in Section 4.2) to increase the S/N and dynamic range 
in UV luminosity at which we can still reliably derive 
typical properties of sources (e.g., stellar mass). In the 
simulation, we find that via stacking, we can achieve 
this 1.5-2.0 mag further in UV magnitude, down to 
Afuv ~ —(18.0-18.5) at z =4-6 (hereafter M t h re sh 2 (•£))■ 
Afthresh 2 (z) corresponds roughly to the UV luminosity of 
the stacked SEDs with S/N ~ 1 at 3.6 /im, below which 
the inferred mass remains uncertain just as one might 
expect. 

In short, our simulation enables us to assess the bias 
and uncertainty of the observed M*-Muv distribution, 
which has so far been generally overlooked in the lit¬ 
erature. We find that the observed M*-Muv distribu¬ 
tion and the inferred relation at high redshift are very 



Fig. 8 .— Probability density function of the recovered M*— Muv 
normalization as a function of intrinsic normalization at 2 : = 7 and 
2 : = 8 in our mock galaxy simulations, when the slope is fixed to 
the intrinsic slope. When the slope is fixed, the normalization can 
be recovered at high confidence. 

sensitive to the choice of stellar mass estimator and the 
UV magnitude range even when using the same data set 
and can be dominated by systematics if not tested thor¬ 
oughly. From this simulation, we derive the redshift- 
dependent UV magnitude thresholds, Mthreshi (z) and 
M t hresh 2 (z), above which we can rely on the median 
mass and scatter of individual galaxies in each UV 
magnitude bin and the median mass of stacks, respec¬ 
tively. Specifically, we find Mthreshi (z) to be Muv = 
—20.0,—20.0,—20.0,—22.5,—22.5 and M t hresh 2 (z) to be 
M uv = -18.0, -18.5, -18.5, -20.5, -20.5 at z = 4, 5, 6, 
7, 8. 

Based on our findings, we derive the best-fit M*-Muv 
relation by combining the median mass of galaxies in 
each UV magnitude bin at Muv < Mthreshi (z) and the 
median mass of stacked points at Mthresbi (z) < Muv < 
M t hresh 2 (z), neglecting galaxies fainter than M t hresh 2 (z) 
when fitting this relation. To explore the validity of this 
optimized method of combining individual bright galax¬ 
ies with stacked faint galaxies, we show in Figure 7 the 
distribution of the recovered slope at z = 4-8 for var¬ 
ious input slopes with our optimized method, demon¬ 
strating that we can recover the intrinsic relation fairly 
well even at z = 6. From this simulation, we derive the 
pdf of intrinsic slope (slopes) of the M*-Muv relation 
given the observed slope (slope ou t), based on the results 
from each of the 50 simulations at each input slope. We 
find, for the given slope observed from our real sam¬ 
ple, the central 68% range of the intrinsic slope to be 
—0.57 < slopein < —0.63, —0.42 < slopes < —0.53, and 
—0.42 < slopejn < —0.64 at z = 4, 5, and 6, respec¬ 
tively. This is comparable with the la conhdence level 
of the observed slope, indicating that our M*-Muv rela¬ 
tions and their uncertainties at z = 4-6 measured from 
our real sample are relatively reliable, when we use this 
optimized method. 

It is encouraging that we find no bias at z = 7-8, but 
the broad pdf of the recovered slope for a given input 
slope at z = 7-8 indicates that it is hard to constrain the 
intrinsic slope with the currently available data. There¬ 
fore, we perform another test to see if our current data 
can constrain the normalization of the intrinsic M*—Muv 
relation when we apply an additional prior of a known 
intrinsic slope. Figure 8 shows the pdf of the recovered 
normalization as a function of intrinsic normalization at 
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TABLE 1 

Best-fit Fiducial log 10 (M*)-M uv relation 


z 

Normalization 

(M e ) 

Slope 

log M*( Muv = _ 2 i) 
(Me) 

4 

-1.70 ± 0.65 

-0.54 ± 0.03 

9.70 ± 0.02 

5 

-0.90 ± 0.74 

-0.50 ± 0.04 

9.59 ± 0.03 

6 

-1.04 ± 0.57 

-0.50 ± 0.03 

9.53 ± 0.02 

7 

-1.20 ± 0.16 

-0.50 ± — 

9.36 ± 0.16 


(-4.23 ± 2.12) 

(-0.65 ± 0.10) 

(9.45 ± 0.07) 

8 

-1.56 ± 0.32 

-0.50 ± — 

9.00 ± 0.32 


(-5.47 ± 6.19) 

(-0.69 ± 0.31) 

(9.00 ± 0.31) 


Note. — At 2 ~ 7—8, numbers in parentheses are the best- 
fit parameters when we do not fix the slope to the best-fit slope 
at 2 ~ 6. The quoted errors represent the la uncertainties. The 
normalization is defined as the logarithmic stellar mass at Muv = 

0 (logM, (Muv=0) )- 

z = 7-8 when the slope is fixed to the intrinsic slope, 
illustrating that the normalization can be recovered ac¬ 
curately within ±0.5. 

The best-fit relation for our real sample in Section 4.4 
is derived with the optimized method of combining in¬ 
dividual bright galaxies with stacks of fainter galaxies. 
When we need to use the full M*-Muv distribution or a 
scatter for our GSMF derivation, we use those inferred 
at Muv < Mthreshl(z)- 

4.4. M*-Muv Relation 

Using our optimized method vetted in Section 4.3, we 
derive the best-fit linear log lo (M*/M 0 )-Muv relation at 
each redshift. Specifically, we use the median mass of in¬ 
dividual galaxies for bins with Muv < M t hreshi(z) and 
the median mass of median flux-stacked SEDs for bins 
with Afuv < A7 t hresh2 ( 2 ) • Our best-fit mass-to-light rela¬ 
tion, listed in Table 1, has a slope of —(0.50-0.69), which 
is slightly steeper than a constant mass-to-light ratio of 
a slope of —0.40. As the slope between 2 = 4 and z = 6 
is nearly constant (~ —0.5) but the uncertainty at 2 — 
7-8 is large, we fix the slope at z = 7 and z = 8 to be 
the same as the best-fit slope at z = 6 while leaving the 
normalization as a free parameter. 18 Our mock galaxy 
simulation in Section 4.3 indicates that our current data 
allow us an accurate recovery of the intrinsic normal¬ 
ization at z = 7-8 with a prior of a known input slope 
(Figure 8). We find a (very) weak evolution in the nor¬ 
malization between z = 4 and z = 7 (a decreasing nor¬ 
malization with increasing redshift). The normalization 
evolves from log(M*/M 0 )(m uv =_ 2 i) = 9-70 (at z =4) to 
9.36 (at z = 7) at only 2cr significance. Interestingly, 
the normalization appears to decrease more rapidly from 
z = 7 to 2 = 8 than at lower redshifts, although the small 
sample at z = 8 prevents drawing any firm conclusions. 
We discuss this further in Section 6. 

Although the M*-Muv distribution of our flux-limited 
sample has nonzero scatter, the derived Af*-Muv rela¬ 
tion is not subject to Malmquist bias (i.e., missing faint 
galaxies at a given stellar mass) as we estimate the Af*- 
Muv relation in bins of luminosity and not stellar mass. 
Therefore, the derived relation should be robust against 
the Malmquist bias, which could artificially result in a 

18 Albeit slightly shallower (—(0.44—0.48)), the slope predicted 
from the SAMs of Somerville et al. (2016, in preparation) is almost 
constant over the redshift interval as well. 


steeper slope than the intrinsic one by losing the faint 
envelope of galaxy distribution for a given stellar mass. 19 

There are discrepancies of 0.3-0.7 dex between differ¬ 
ent studies in the measured median mass at a given UV 
magnitude even at z ~ 4, which are larger at fainter 
UV bins (Gonzalez et al. 2011; Lee et al. 2012; Stark et 
al. 2013; Duncan et al. 2014; Salmon et al. 2015). This 
may reflect a number of systematic uncertainties associ¬ 
ated with sample selection and stellar mass estimation. 
As these uncertainties make a direct comparison between 
different studies difficult, it highlights the importance of 
a comprehensive and independent analysis to verify sys- 
tematics. 

Overall, we found the best-fit slope shallower than that 
of Gonzalez et al. (2011) but similar to those of Stark 
et al. (2013) and Duncan et al. (2014), with an excep¬ 
tion at 2 = 4. First, the slope of our best-fit relation 
of —0.54(±0.03) at 2 = 4 is significantly shallower than 
that of Gonzalez et al. of —0.68(±0.08), which is based 
on an order-of-magnitude-smaller sample. As the rela¬ 
tion of Gonzalez et al. is derived with no nebular cor¬ 
rection, it is not surprising that their stellar masses for 
UV-bright galaxies are higher than ours. However, their 
stellar masses for galaxies in faint UV bins are lower than 
ours by ~0.2 dex at Muv ~ —18, resulting in a steeper 
slope than ours. Meanwhile, the M*-Muv relation of 
Stark et al. (2013) shows a good agreement with ours 
at all redshifts with only a slightly shallower slope than 
ours. Duncan et al. (2014) in general find higher stel¬ 
lar masses for galaxies than we do in faint UV bins, re¬ 
sulting in a higher normalization and a shallower slope 
than ours, with the biggest difference in normalization 
(Alog(M*/M 0 )M uv =- 2 i = 0.25) and slope (A = 0.09; 
2.5cr significance) being observed at 2 = 4 (see Section 
5.4 for more discussion). 

5. THE GALAXY STELLAR MASS FUNCTION 
5.1. Derivation of the GSMFs 

We now convolve our M*-Muv distribution with the 
observed rest-frame UV luminosity function to derive 
GSMFs. The luminosity functions we utilize in this 
study are from Finkelstein et al. (2015), which included 
the full CANDELS GOODS fields, the HUDF and two 
parallel fields, and the Abell 2744 and MACS J0416.1- 
2403 parallel fields from the Hubble Frontier Fields data 
set. These luminosity functions are already corrected for 
incompleteness and selection effects using the detection 
probability kernels derived from fake source simulations, 
as discussed by Finkelstein et al. (2015). 

Figure 9 presents GSMFs constructed using four dif¬ 
ferent methods: 

1. “Raw Bootstrapped GSMF”: 

We constructed the “observed” GSMF by combining 
the UV luminosity function with the observed M*-Muv 

19 As discussed in Section 4.1, objects with high mass and low 
UV flux that are observed or that we may be missing are believed to 
be the subdominant population for all UV luminosity bins and thus 
would not change the derived relation (see also Salmon et al. 2015 
for a similar argument based on the tight scatter of the star-forming 
main sequence). When probing in bins of stellar mass, however, 
those massive outliers seen at z = 4-5 with log> 10 are 
off from the linear correlation derived in bins of UV luminosity in 
this section. 
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Fig. 9.— Our fiducial galaxy stellar mass functions, from a J + H -band selected sample of galaxies at 2 : = 4—8 (red filled circles). The 
open squares, filled squares, small stars, and small red filled circles indicate raw bootstrapped, incompleteness-corrected bootstrapped, 
constant-scatter, and asymmetric-scatter GSMFs, respectively. The red dotted, red dashed, red dash-dot and red solid lines represent the 
Schechter fit for the last three and for our fiducial GSMFs, respectively. The uncertainty on our fiducial GSMF includes contributions 
from both the UV luminosity function uncertainties and the uncertainty in the relation. The light-red shaded regions denote 

1,000 Schechter fits for our fiducial GSMF randomly chosen within the lcr 3-dimensional contour of the Schechter parameters determined 
from our MCMC analysis. The blue points (open squares, filled squares, dashed line) correspond to previous estimates (raw bootstrapped, 
incompleteness-corrected bootstrapped, constant-scatter GSMF) of Gonzalez et al. (2011) from WFC3/IR data of the ERS (for 2 : = 4-6) 
and from WRC3/IR data of the ERS, HUDF09, and NICMOS data over the GOODS fields (for z = 7). Also overplotted are recent 
estimates for GSMFs from the literature—from As-band selected sample at 2 : ~ 4 (Ilbert et al. 2013; Muzzin et al. 2013), from 4.5/j.m 
selected sample (Caputi et al. 2011, 2015), and from rest-frame UV-selected samples (Stark et al. 2009; Santini et al. 2012; Lee et al. 2012; 
Duncan et al. 2014; Grazian et al. 2015). All points and lines are converted to a Salpeter IMF. The thick grey lines show dark matter 
halo mass functions scaled to a baryon conversion efficiency of 20% (i.e., 20% of halo mass times the cosmic baryon fraction of 
Our GSMFs are characterized by a steeper low-mass-end slope of —1.55 (—2.25) at 2 : = 4 (z = 8) compared to that of Gonzalez et al. of 
-1.43 (-1.55) at z = 4 (z = 7). 
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distribution. We first drew 10 5 galaxies from the best- 
fit Schechter function of the UV luminosity function in 
the range of —30 < Myv < —13. Then, we assigned 
stellar mass for each galaxy to be the stellar mass of 
the randomly chosen galaxy with similar rest-frame UV 
absolute magnitude from the observed A/*-A/uv distri¬ 
bution. This method accounts for outliers, which may 
be a non-negligible fraction of galaxies at z = 4. 

The non-negligible spread in stellar mass at a fixed 
rest-frame UV absolute magnitude observed in Figure 5 
can result in GSMFs with an underestimated low-niass- 
end slope if we do not account for the “unobserved” pop¬ 
ulation of galaxies below the detection threshold of the 
survey. This is because galaxies at a given UV luminos¬ 
ity have a range of mass-to-light ratios; thus, galaxies be¬ 
low the detection limit can still have stellar masses high 
enough to contribute to the number density in the stellar 
mass bins of our interest (affecting the last few points in 
the GSMFs depending on the amount of spread). There¬ 
fore, we need to correct for the incompleteness by as¬ 
suming an A/*-A/uv distribution in UV luminosity bins 
below the current sensitivity limit. We use three different 
incompleteness correction schemes. 

2. “Incompleteness-corrected Bootstrapped GSMF”: 

A reasonable assumption on the A-:/*-A/uv distribution 

in the unobserved faint UV luminosity bins is that the 
observed distribution of individual galaxies in bright UV 
bins above a threshold represents the intrinsic distribu¬ 
tion and can be extended to fainter bins. In Section 4.3, 
we found the redshift-dependent threshold M t hreshi (z) at 
each redshift above which the intrinsic distribution can 
be recovered well without any noticeable bias or increase 
in scatter. We extend the observed A/*-Aiuv distribu¬ 
tion at Afuv < A--fthreshi(^) to fainter UV luminosities 
down to A/uv = —13, keeping the distribution centered 
around the best-fit A/*-Muv relation. The rest of the 
procedures are the same as those for the “raw boot¬ 
strapped GSMF”. 

3. “Constant-scatter GSMF”: 

Instead of using individual points in the AT*-A/uv dis¬ 
tribution, we assumed an idealized lognormal distribu¬ 
tion around the best-fit A/*—A/uv relation with a con¬ 
stant scatter, inferred from bright UV luminosity bins 
that have a statistical number of galaxies and high com¬ 
pleteness. We used the mean scatter estimated from the 
two faintest bins at Muv < Afthreshi (z) with more than 
10 galaxies, which is ~0.4-0.5 dex. The UV luminosity 
function was then convolved with this lognormal distri¬ 
bution to derive the “constant-scatter GSMF”. 

4. “Asymmetric-scatter GSMF”: 

The three approaches above are basically the same as 
those used by Gonzalez et al. (2011). However, the mass 
distribution of our sample in a given rest-frame UV abso¬ 
lute magnitude bin is not symmetric with respect to the 
best-fit relation. Rather, the lower side of the best-fit re¬ 
lation has in general a smaller scatter. As already noted 
in Section 4.1, the lack of galaxies with high UV luminos¬ 
ity and low mass contributes in part to the asymmetric 
scatter, and the results in Section 4.3 in addition indi¬ 


cate that it could be an intrinsic property and not just 
an observational bias. Therefore, we assume a lognormal 
distribution with an asymmetric scatter with respect to 
the best-fit A/*-A/uv relation (a different sigma above 
and below the mean) inferred from the two faintest bins 
at A/uv < Afthreshi (z) and extend it to fainter UV lu¬ 
minosity bins. 20 The resultant GSMF constructed via 
this method is referred to as the “asymmetric-scatter 
GSMF”. 

While the asymmetric-scatter GSMFs were devised to 
probe the incompleteness-corrected low-mass-end slope, 
they might not account properly for the fraction of out¬ 
liers seen at z = 4 and z = 5 that can impact the 
high-mass end of the GSMF. Thus, we combine the 
bootstrapped GSMF for high masses (log(A/*/A/ 0 ) > 
10) and the asymmetric-scatter GSMFs for low masses 
(log(Af* /Mq) > 10), and consider this our fiducial 
GSMF. We found that the GSMFs computed using 
the 1/Umax method (Schmidt 1968, gray diamonds in 
Figure 9) are in excellent agreement with our fidu¬ 
cial GSMFs at the high-mass end. Figure 9 shows 
our fiducial GSMFs at z = 4-8 (listed in Table 2) as 
well as those derived from the four methods described 
above. The incompleteness-corrected, constant-scatter, 
and asymmetric-scatter methods yield consistent results 
(except at z = 8 where the GSMF is less robust due 
to the small sample size), and their Schechter fits (see 
Section 5.2) are nearly indistinguishable within the un¬ 
certainties on the low-mass end. 

5.2. Schechter Fit and Uncertainties 

Before we parameterize our GSMF with a Schechter 
function, we first need to derive the uncertainties for our 
GSMF data points. Randomly perturbing points of the 
GSMFs is not a proper way to estimate the uncertainties 
of the GSMFs because they are correlated. Thus, we de¬ 
rive the 68% confidence interval of our fiducial GSMFs 
as follows. First, we randomly drew 1000 samples from 
a Markov Chain Monte Carlo (MCMC) chain of the 
UV luminosity function Schechter parameters derived by 
Finkelstein et al. (2015) within the 3-dimensional lcr con¬ 
tour of (L*,aL,(f>* L ). Then, for each luminosity function 
generated from the Schechter parameters, we assigned an 
A/*-Muv relation with (slope, normalization) values ran¬ 
domly chosen within the 2-dimensional lcr contour of the 
best-fit parameters of the Af*-Afuv relation. 21 The new 
A/*-Muv distribution was then combined with the new 
luminosity function to generate a GSMF in the same way 
that our GSMF was constructed. This generates 1000 
GSMFs, of which the minimum and maximum represent 
the lcr upper and lower limits of the GSMFs, respectively. 

We parameterize our GSMFs with a Schechter (1976) 
function, 

</>{M,)dM. = (4>*/M*)x(M*/M*) a exp[—(A/* /M*)\dM* 

( 1 ) 

20 For z = 7 and z = 8 where the scatter is not robustly mea¬ 
sured due to the small sample size, we assume that the 
distribution follows that at 2 : = 6. 

21 For z =7-8 where we fix the slope of the M*-Muv relation 
to the z = 6 value, the uncertainty in the normalization is taken 
into account. 
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which is characterized by a power law with a low-mass- 
end slope of a, an exponential cutoff at stellar masses 
larger than a characteristic mass, M*, and a normal¬ 
ization (j)*. The best fit, uncertainties, and posterior 
distributions of the Schechter parameters for our fidu¬ 
cial asymmetric-scatter GSMFs were derived by running 
an MCMC algorithm that samples the 3-dimensional pa¬ 
rameter space of the Schechter parameters. To ensure 
full coverage of the parameter space and assess conver¬ 
gence, we ran five parallel chains composed of 10 5 steps 
each. The starting position of the chain was determined 
by a coarse grid search that minimized the x 2 statis¬ 
tic between the model and the data. The first 10% of 
steps were disregarded in the burn-in phase before run¬ 
ning each chain to reduce the dependence of the poste¬ 
rior distribution on the initial position. The proposal 
distribution of each parameter was assigned as a nor¬ 
mal (lognormal) distribution for a (M* and (f>*) with 
the width tuned to generate an acceptance rate of ~ 
23-30%. 22 As a prior, we limited the sampling param¬ 
eter space to be a > —10, 8 < log (M*/Mq) < 13, 
and log(0*/Mpc -3 ) > —8. For 6 < z < 8, where the 
constraints on M* are weak (see the large error bars 
on the open gray circles in Figure 10), we took a log¬ 
normal prior on M* with mean log(M*/M 0 ) = 10.75 
and standard deviation of 0.3 dex, following the poste¬ 
rior distribution of M* at z = 4 and z = 5. Comparing 
the likelihood of the current step, which is defined as 
£ oc e~ x / 2 , with the proposed set of parameters deter¬ 
mines whether the proposal is accepted. We employ the 
Metropolis-Hastings algorithm for the acceptance crite¬ 
ria. After running all chains, convergence was assessed by 
examining trace plots of parameters, as well as using the 
Rubin-Gelman R diagnostic (Gelman & Rubin 1992) for 
each marginalized posterior distribution. The diagnos¬ 
tic value R ~ 1 suggests convergence, and we confirmed 
that for all redshifts and parameters the diagnostic has 
a value 1.00 < R < 1.01. 

From the resulting joint posterior distribution, we ex¬ 
tracted the marginal posterior distribution of each pa¬ 
rameters. The median and the central 68% of the 
marginal posterior distributions provide our fiducial 
Schechter parameters and an estimate of the 68% confi¬ 
dence interval on each parameter (shown as the error bars 
in Figure 10). In short, the uncertainties on the GSMFs 
and Schechter parameters include (i) the uncertainty of 
the best-fit M*-Muv relation and (ii) the uncertainty of 
the Schechter parameters in the UV luminosity function, 
the latter of which includes Poisson errors. Other sources 
of random errors on the derived GSMFs, including cos¬ 
mic variance, are discussed in Section 7.2. 

Our best-fit Schechter parameters are listed in Table 
3 and plotted in Figure 10 as a function of redshift. 
Our data support a decreasing (steepening) of the low- 
mass-end slope (a ~ -1.55±g;[$, -1.70^;^, -1.9lt^9, 
— 1.951 q}®, — 2 . 25+035 at z = 4, 5,6,7, 8) with increas¬ 
ing redshift, asymptoting to the faint-end slope of the 
UV luminosity function (a^ ~ —1.56ig'og, — l-67jj{}'o6> 

22 The ideal acceptance rate, which affects the efficiency of the 
chain, depends on the dimension ( d ) of the problem, decreasing 
from ~44% for d = 1 to ~23% for d oo (Gelman et al. 1996, pp. 
599-607). 



Fig. 10. — Redshift evolution of the best-fit Schechter parameters 
for our fiducial GSMFs. The low-mass-end slope, a, evolves toward 
a steeper value with increasing redshift, asymptoting to the faint- 
end slope of the UV luminosity function. Conversely, the open 
gray circles denoting the best-fit Schechter parameters with a flat 
prior on M* show that we observe no significant evolution in the 
characteristic mass, M*, though our observations do not allow us 
to place reasonable constraints on M* at z > 6. 



Fig. 11. — Confidence contours of the best-fit Schechter parame¬ 
ters for our fiducial GSMFs at z = 4-8 at the 68% and 95% levels, 
showing that we can place reasonable constraints on a. The best- 
fit Schechter parameters at z = 6—8 were derived with a lognormal 
prior on M*, and their contours are denoted as dot-dashed lines. 
The best-fit values are shown as filled circles. 
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-2.02t°;}°, M2.03lo.2O’ - 2 - 36 -o:40 at z = 4,5,6, 7,8), 
and possibly a decrease in </>* with increasing redshift as 
well. Conversely, our GSMFs favor no evolution in Al* 
with redshift, although we cannot rule out the possiblity 
of evolution due to the large uncertainties, in particular 
at z > 6 (see the gray open circles in Figure 10, which 
are derived with a flat prior on Al*). The large error 
bars on Al* are because constraints on the bright end of 
the luminosity function (and thus the massive end of the 
GSMF) remain weak, as even our wide-area data reach 
only 1-2 mag brighter than the characteristic magnitude 
of the UV luminosity function, L* (marked as arrows at 
the bottom of Figure 5), and the UV magnitude bins 
brighter than L* are populated by only a few galaxies. 
The lack of robust constraints on the massive end of the 
GSMF leads to a well-known degeneracy between Al*, a , 
and <j>*, as shown in Figure 11 with confidence contours 
on the Schechter parameters. 

5.3. Eddington Bias 

Several recent studies (Ilbert et al. 2013; Caputi et 
al. 2015; Grazian et al. 2015) investigated the effect of 
the Eddington bias (Eddington 1913) on the low-mass- 
end slope of the GSMF, finding that the increase in the 
stellar mass uncertainties for individual galaxies at low 
masses results in an artificial steepening of the low-mass- 
end slope (Grazian et al. 2015; Caputi et al. 2015). In our 
methodology, the “effective” uncertainty in stellar mass 
is not stellar mass dependent, as we do not use individual 
stellar masses to construct the GSMF on the low-mass 
end but rather convolve the UV LF with stellar mass-UV 
luminosity distribution, where we use the distribution 
around ~L* (—21 < Afuv < —20) to extend to fainter 
luminosities. Taking a similar approach to other studies, 
we investigated the effect of the Eddington bias by pa¬ 
rameterizing the observed GSMF as a convolution of an 
“intrinsic” Schechter function with a lognormal function 
with its width to be the mean stellar mass uncertainty 
of galaxies with —21 < Muv < —20. We implemented 
a grid search for the “intrinsic Schechter parameter to 
correct for the Eddington bias, finding that while we see 
a similar trend to what Grazian et al. (2015) found (the 
“intrinsic” low-mass-end slope is similar [at z = 4-5] or 
slightly flatter [at z > 6], driven by the change in Al*, 
which is uncertain), the change in the number density on 
the low-mass end is negligible and the trend of steepen¬ 
ing low-mass-end slope with increasing redshift persists 
after correcting for the Eddington bias. 

5.4. Comparison with Previous Work 

Figure 9 compares our GSMFs with recent estimates 
from the literature determined by rest-frame UV-selected 
galaxies (either using photometric redshift or color-color 
selection; Stark et al. 2009; Gonzalez et al. 2011; Lee et 
al. 2012; Santini et al. 2012; Duncan et al. 2014; Grazian 
et al. 2015) 23 or rest-frame optical-selected galaxies (Ca¬ 
puti et al. 2011; Ilbert et al. 2013; Muzzin et al. 2013; Ca¬ 
puti et al. 2015). All GSMFs are converted to a Salpeter 
IMF when necessary. Overall, the comparison of our 
GSMFs with previous estimates demonstrates that there 

23 For the GSMFs of Stark et al. (2009), we apply correction 
factors for nebular emission of xl.l, 1.3, and 1.6 at z = 4, 5, 6 that 
are inferred in Stark et al. (2013). 


exists a considerable discrepancy in the high- and/or low- 
mass ends and in the normalization of the GSMFs at all 
redshifts. However, before discussing the discrepancy in 
detail, we point out that the error bars of the GSMFs 
compared in Figure 9 include only random uncertain¬ 
ties. For example, those of Grazian et al. (2015) include 
Poisson errors and errors (photometric scatter and pho¬ 
tometric redshifts) derived from their Monte Carlo simu¬ 
lations, while those of Duncan et al. (2014) include Pois¬ 
son errors and photometric redshift uncertainty. Those 
of Gonzalez et al. (2011) include the uncertainties of the 
UV luminosity function of Bouwens et al. (2007, 2011) 
and the scatter of the Af*-Afuv relation. None of the 
plotted points include cosmic variance or other sources 
of systematic uncertainties. 

A major source of discrepancy in the GSMFs between 
different studies may be attributed to the systematic un¬ 
certainties associated with stellar mass estimation, as 
already noted by many studies (e.g., Marchesini et al. 
2009; Mobasher et al. 2015). Each study adopts different 
assumptions on the SFH, dust law, metallicity, nebular 
emission, etc., in their SPS modeling, of which effects 
on the derived stellar mass can be significant. More¬ 
over, many parameters are degenerate, thus making as¬ 
sessment of the systematic effects induced by these dif¬ 
ferent assumptions on the disagreement of the observed 
GSMFs hard to achieve. While investigating the system¬ 
atic effects caused by the different sets of assumptions 
adopted in previous studies is beyond the scope of our 
study, we stress that we have focused on deriving the in¬ 
trinsic GSMFs by exploring the systematic effects inher¬ 
ent in our analysis, minimizing them via our SAM+mock 
galaxy simulations (see Section 4.3). 

With these above caveats in mind, we now discuss 
the discrepancy highlighted by the direct comparison be¬ 
tween the GSMFs in Figure 9. First, we observe a dis- 
greement in the normalization of the GSMFs between 
ours and the estimates from the literature at all red¬ 
shifts. At z = 4 and z = 5, although our normalization 
is in good agreement with Lee et al. (2012) and Grazian 
et al. (2015), a prominent discrepancy is observed with 
Gonzalez et al. (2011), Santini et al. (2012), and Duncan 
et al. (2014). The GSMFs of Gonzalez et al. (2011) are 
found to lie systematically below ours and others in the 
literature at z = 4 and z = 5, and the GSMFs of Santini 
et al. (2012) at z ~ 4 and Duncan et al. (2014) at all red¬ 
shifts are found to lie above. At higher redshifts of z = 
6-7, the normalization of our GSMFs is lower with re¬ 
spect to those of Duncan et al. (2014) and Grazian et al. 
(2015, at z = 6) but shows a good agreement with Stark 
et al. (2009) (at the massive end; they do not correct for 
incompleteness), Gonzalez et al. (2011), and Grazian et 
al. (2015, at z = 7). Although the overall difference in 
the normalization of the GSMFs between different stud¬ 
ies remains similar at z = 7, the larger error bars at z = 7 
render any differences at that high redshift insignificant. 

Examining the different mass regimes, we notice an in¬ 
teresting disagreement at z = 4, where wide-area ground- 
based surveys (Ilbert et al. 2013; Muzzin et al. 2013), 
which are potentially more sensitive to more rare, mas¬ 
sive galaxies, may be more robust. Specifically, in the 
most massive bin of our study at log(Af*/A Iq) = 11.25, 
both Muzzin et al. (2013) and Ilbert et al. (2013) found 
a number density ~ 0.1-0.3 dex higher than we find 
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TABLE 2 

Galaxy stellar mass function atz= 4-8 


log M* 



logc/> 



(Me) 



(dex -1 Mpc -3 ) 




z = 4 

z = 5 

z = 6 

z = 7 

z = 8 

7.25 

i C7+O.2I 
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1,41 ' -0.32 

— 1 03+ 0 - 54 
l.OO —0 54 

I 7q+1.01 
J -- < °-0.84 

7.75 

1 77+0.15 

1 '-0.14 

1 79+O.2O 
-‘-•'^-0.20 

1 81 +0-23 
i -° i -0.28 

—2 07 +0-45 
z - u '-0.41 

9 90+O.84 
^.zo_o.64 
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—2 QO +0 ' 13 
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Z ' U± -0.16 
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z ’^ y -0.32 
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8.75 
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9 oq+0.15 

z.oo-o 10 
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Z - Ot> -0.15 

—2 96 +0 ‘ 32 
z,yD -0.30 

—3 45+ 0 - 57 
°’ q:C> -0.60 

9.25 
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^•°^_0.09 

-2-68 ±°;°I 

— 3 14 +0-12 

0.14_o 11 

— 3 47+ 0 - 32 
-0.35 

—4 21 +°- 63 

4 ’ Z1 -0.78 

9.75 

—2 qi+ 012 
z ' yi -0.05 

q -1 9 + 0.09 
O.-L^_0.11 

— 3 69 +0 12 
° lUy —0.13 

—4 ii+ 0 - 41 

4 - ±J --0.57 

— 3 31 + 101 

O.Ol-i 64 

10.25 

o 07+O.O9 
0,0 ' -0.12 

—3 47+ 0 - 16 
-0.14 

—4 27 +0 ' 38 
4 - z '-0.86 

—4 6i+ 0 - 72 
4 - D1 -0.82 

— 

10.75 

—4 QO +0 ' 20 

4 ' uu -0.25 

-4.12±°-ji 

— 

— 3 94+ 0 - 90 
D - Z4 -0.57 

— 

11.25 

—4 34+0-34 
4.04 —0 55 

-4.88i°;« 

— 

— 
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Note. — The quoted la errors include the uncertainties of the UV luminosity 
function and the M*-Muv relation. 


TABLE 3 

Best-fit Schechter function parameters of 
OUR FIDUCIAL GSMFS 


z 

log M* 

a. 

4> * 


(Me) 


(10“ 5 Mpc" 3 ) 

4 

10.50±°; 3 ° 

1 33+O.O8 
± -°0 —0.07 

23 fi8"*" 21 ' 75 
zO.Oo_ 12 g0 

5 

io.97±°;“ 

— 1 70+°-° 8 
± -' u -0.07 

3 16~*~ 7 ' 05 
°- iO -4.08 

6 

10.72«; 2 ® 

— 1 ai+°- 09 
0.09 

1 qc + 1.66 
l.OO —0.75 


(11.08 j+Q4) 

(-l.90±g;lg) 

(^ + - 0 H) 

7 

10.78i°;l 

-1.95±g;ii 



(11.33±l°f) 

(-l.getg; 2 ?) 

(0.l3±g;f^) 

8 

10 72+ 0,29 
1U. 1 z — 0 29 

9 93+0-72 
Z - Z0 —0.35 

0-035±g;gSg 


(10.301); 31 ) 

(—2-l2±g;|g) 

(0.125™) 


Note. — The quoted best-fit values and la errors of 
the Schechter parameters represent the median and the 
central 68% confidence interval of the marginal posterior 
distribution of each parameter obtained from our MCMC 
analysis. At 6 < z < 8, we show in parentheses the 
results derived with a flat prior on M*. 

(and ~0.2-0.3 dex higher than Duncan et al. (2014) and 
Grazian et al. (2015)). This discrepancy may be at¬ 
tributed in part to the fact that the median redshift of 
these ground-based surveys is z = 3.5, lower than ours, 
as well as the other studies shown (z ~ 4). Moreover, 
the former is derived from a shallow but wide ( 1.6 deg 2 
down to K s ~ 24) K s -band (rest-frame optical) selected 
catalog. Thus, they should be more complete and less 
susceptible to cosmic variance and Poisson noise at the 
high-mass end than other works that are based on small- 
field, rest-frame UV-selected catalogs. Meanwhile, Ca- 
puti et al. (2015), who investigated the high-mass end of 
the GSMFs by adding the contribution of -ft+band faint 
but 4.5/im bright ([4.5] < 23) galaxies to the previous de¬ 
terminations of the /F s -band-selected GSMFs by Caputi 
et al. (2011) and Ilbert et al. (2013), present consistent 
results with ours, which may be attributed to the lower 
normalization of the GSMFs of Caputi et al. (2011) with 
respect to others. 

At the high-mass end, GSMFs based on rest-frame UV- 
selected galaxies (Stark et al. 2009; Gonzalez et al. 2011; 


Lee et al. 2012; Duncan et al. 2014; Grazian et al. 2015) 
agree reasonably well with each other when the cosmic 
variance is accounted for (see Section 7.2.2). However, 
the high-mass end at z = 4 from Santini et al. (2012) 
still shows a mild tension with ours due to their overall 
higher normalization with respect to others. 

Turning to the low-mass end, while our survey volume 
is smaller than those of the ground-based surveys of Il¬ 
bert et al. (2013) and Muzzin et al. (2013), the deep data 
set in this study enables us to reach lower in mass than 
these surveys can (log(M*/Mg) > 10), allowing more ro¬ 
bust constraints on the low-mass-end slope. Starting at 
z = 4, our results at the low-mass end are consistent with 
most previous rest-frame UV-selected studies (Gonzalez 
et al. 2011; Lee et al. 2012; Duncan et al. 2014; Grazian 
et al. 2015), with the exception of Duncan et al. (2014). 

Because Grazian et al. (2015) restrict their analy¬ 
sis to higher masses, the only points for comparison 
at log (M*/Mq) < 9 are those from Gonzalez et al. 
(2011) and Duncan et al. (2014). Interestingly, the 
largest disagreement at the low-mass end is found at the 
lowest redshift of 2 = 4, where Duncan et al. (2014) 
find significantly higher number densities (~0.5 dex at 
log (M„/M 0 ) rsj 9) and a steeper low-mass-end slope 
(a ~ —1.9) with respect to the others. This may re¬ 
sult from differences in the measured M*-Muv relation, 
as well as in the faint-end slope of the UV luminosity 
function between our study and that of Duncan et al. 
(2014). PA z = 4, Duncan et al. (2014) found a shal¬ 
lower M uv slope than what we find here. One differ¬ 
ence in methods is that rather than using our hybrid ap¬ 
proach of using individual high-mass galaxies and stacks 
of lower-mass galaxies, Duncan et al. (2014) fit their 
Af*-Aiuv distribution over a wide stellar mass range 
down to \og(M*/Mo) ~ 8. As shown in their simula¬ 
tions (see their Fig. 5), stellar masses for galaxies with 
\og{M*/M q) < 9 are biased toward higher masses (a sim¬ 
ilar result to what we find here), leading to their deriva¬ 
tion of a shallower M*-ALuv slope. A shallower slope of 
the Af*-Muv relation translates into a steeper low-mass- 
end slope and a higher normalization of the GSMF: for 
a given number density of galaxies in bins of UV lumi¬ 
nosity (<+,), the number density of galaxies in bins of 
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stellar mass (4>m) is given by <j>M oc 4>L(dL/dM), and 
is thus higher for a shallower M*-Muv slope. We do 
note, however, that Duncan et al. (2014) did not use this 
M*-Muv relation to measure their GSMF; they used a 
1/Knax method. However, the biases inherent in mea¬ 
suring masses from individual poorly detected galaxies 
may still be responsible for their steeper low-mass-end 
slope at z = 4. Moreover, their faint-end slope of the 
UV luminosity function at 2 ~ 4 is steeper (by ~0.2) 
than that from Finkelstein et al. (2015), on which our 
GSMF is based. 

Unsurprisingly, differences are thus found in the evo¬ 
lution of the low-mass-end slope with redshift. The 
Schechter fit for our GSMF indicates steeper low-mass- 
end slopes. While Gonzalez et al. found a tentative 
steepening in the low-mass-end slope with increasing red- 
shift, the steepening is mild, from —1.43 at z = 4 to 
— 1.55 at 2 = 7. This is a combined effect of a steeper 
faint-end slope of the updated UV luminosity function 
by Finkelstein et al. (2015) used in our analysis and a 
shallower M*-Muv relation found in this study (Section 
4.2). Grazian et al. (2015) and Duncan et al. (2014) both 
find a nearly constant low-mass-end slope of a = —1.6 
and a = —1.9, respectively, and no evidence of the steep¬ 
ening that we observe. 

At lower redshifts of 2 < 4, the consensus is that the 
characteristic mass does not change but the normaliza¬ 
tion evolves (e.g., Marchesini et al. 2009), although the 
evolution of the low-mass-end slope remains controversial 
as some find no evolution (Marchesini et al. 2009) while 
others find a steepening low-mass-end slope with increas¬ 
ing redshift (for a single Shechter function fit; Kajisawa 
et al. 2009; Santini et al. 2012; Ilbert et al. 2013; Tom- 
czak et al. 2014). In this study, at 2 > 4, the observed 
evolution of the GSMF shows a low-mass-end slope that 
steepens with redshift. Additionally, our results tenta¬ 
tively confirm a roughly constant M *, with a decreas¬ 
ing (j>* with increasing redshift, qualitatively similar to 
results at lower redshift (e.g., Ilbert et al. 2013), though 
we acknowledge that our relatively small volume prevents 
robust constraints on M*. 

6. STELLAR MASS DENSITY 

To measure the stellar mass density at 2 = 4-8, we 
integrated the best-fit Schechter function at each red¬ 
shift from 8 < log(M*/M 0 ) < 13, an often adopted 
interval for stellar mass density estimates at high red¬ 
shift in the literature. Table 4 lists our estimates of the 
stellar mass density along with their la uncertainties, 
calculated as the minimum and maximum stellar mass 
densities allowed by the 3-dimensional la contour of the 
Schechter parameters obtained in Section 5.2. Figure 
12 presents the evolution of the stellar mass density at 
2 = 4-8, alongside values compiled from the literature 
(converted to a Salpeter IMF when necessary). Most 
data points from the literature are taken from the com¬ 
pilation by Madau & Dickinson (2014), with the excep¬ 
tion of Gonzalez et al. (2011), which Madau & Dickinson 
(2014) corrected for nebular emission. We instead show 
the uncorrected (original) points together with Stark et 
al. (2013), which are the Gonzalez et al. values corrected 
for nebular emission. We also add the recently published 
works of Duncan et al. (2014) and Grazian et al. (2015). 
The error bars from most of the studies include only ran- 
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Fig. 12.— Evolution of the stellar mass density. The stellar 
mass densities were obtained by integrating the best-fit Schechter 
functions for our fiducial GSMFs between M* = 10 8 Mq and 
M * = 10 13 Mq (red circles). The error bars indicate the minimum 
and maximum values of stellar mass density allowed by the lcr 
contour of the Schechter parameters. Small symbols are the com¬ 
pilation of the stellar mass densities from the literature by Madau 
&; Dickinson (2014) (using their colors and symbols) along with re¬ 
cent estimates from Stark et al. (2013), Duncan et al. (2014), and 
Grazian et al. (2015), listed in the legend. All points and lines are 
converted to a Salpeter IMF. The solid black curve marks param¬ 
eterization of the time integral of SFRD from Madau Sz Dickinson 
(2014) after gas recycling (R = 0.27) is accounted for, represent¬ 
ing the prediction for the stellar mass density. The stellar mass 
densities predicted for the two scenarios suggested by Oesch et al. 
(2014) for the SFRD evolution at z > 8 are denoted as the black 
and blue dashed lines (see Section 6 for more details). Our stellar 
mass densities show a remarkable agreement with estimates of the 
stellar mass density from the SFRD. Other colored solid curves are 
stellar mass densities predicted from the three SAMs introduced in 
Section 7.1.1, and the cyan hatched region denotes the 95% poste¬ 
rior range of the Lu et al. (2014) model. For reference, we denote 
fractions of the local stellar mass density measurement (Baldry et 
al. 2012) as horizontal dotted lines. 

dom errors. 

Our estimates of the stellar mass density at z = 4-5 are 
in broad agreement with previous measurements within 
the uncertainty, with the exception of Caputi et al. (2011) 
at z = 5 and Duncan et al. (2014) at z = 4. Duncan et al. 
found a ~0.5 dex higher stellar mass density at z = 4, a 
deviation at 1.9cr, mainly due to their steeper low-mass- 
end slope compared to ours. The stellar mass density at 
z ~ 5 from Caputi et al. shows a value lower by about 
0.8 dex than our estimates (a deviation at 4.2cr), which 
is surprising, given that their median redshift is slightly 
lower (z ~ 4.6). However, the sample selection of Caputi 
et al. is very different from ours, such that their sample 
is selected in the IRAC 4.5 /im band, which, while more 
complete for very red galaxies, may underestimate the 
incompleteness to star-forming galaxies. 

At higher redshifts of z = 6-7, our measurements are 
~0.5-0.7 dex lower than those of Gonzalez et al. (2011) 
and Duncan et al. (2014) but show an excellent consis- 
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TABLE 4 

Cosmic Stellar Mass Density at z = 
4-8 


z 

log p« 


(M® Mpc -3 ) 


7.17 

+0.22 

-0.13 

O 

Q_ 

7.09 

+0.28 

-0.19 


6.48 

+0.31 

-0.18 

X 

CD 

(6.53 

+0.55 \ 

—0.26' 

~o 

6.19 

+0.62 

-0.40 


(6.24 

+0.88 \ 

-0.45' 


5.50 

+0.83 

-0.81 


(5.51 

+1.64\ 

-1.80/ 



Note. — Stellar mass estimates by inte¬ 
grating the best-fit Schechter functions for 
our GSMFs over 8 < log(M*/Mg) < 13. 
The quoted lcr error bars represent the min- 
inum and maximum values of stellar mass 
density allowed by the 3-dimensional 1<7 con¬ 
tour of the Schechter parameters. The values 
in parentheses at 6 < z < 8 are the results 
from the best-fit Schechter function derived 
with a fiat prior on M*, displaying negligi¬ 
ble difference from our fiducial stellar mass 
estimates. 


tency with the recent estimate of Grazian et al. (2015). 
The agreement in the stellar mass density in spite of the 
difference in the normalization of the GSMFs between 
our study and Grazian et al. (2015) at z = 6 (their higher 
normalization compared to ours) is a consequence of their 
shallower low-mass-end slope by 0.4, which compensates 
the difference in the stellar mass density due to their 
higher normalization when integrating the GSMFs. As 
the GSMFs are less well constrained at these high red- 
shifts, our stellar mass density is only in mild tension 
with other studies at z = 6 (< 2cr). At z = 7, although 
the estimates of the stellar mass density from different 
studies differ up to 0.7 dex, the increased error bars mean 
that these differences are not currently statistically sig¬ 
nificant. 

As the stellar mass is to first order the time integral of 
past star formation activity, a comparison of the stellar 
mass density with the time integral of the SFRD should 
yield similar values if both estimates are accurate. At 
high redshifts, however, both quantities have large uncer¬ 
tainties. The limiting factor in determining an accurate 
SFRD is a determination of dust attenuation, for which 
the observed UV luminosity density is corrected, without 
a direct observation of the dust-obscured star formation 
for most cases. Likewise, the uncertainties involved in 
the determination of the stellar mass density (the sys- 
tematics in stellar mass estimates, the uncertainties on 
the abundance of low-mass galaxies, etc.) can impact the 
stellar mass density and potentially result in a mismatch 
between the integral of the SFRD and the stellar mass 
density. 

Figure 12 compares the time integral of the SFRDs 
with the stellar mass density derived in this work. In 
Figure 12, the black solid line indicates the parameter¬ 
ization by Madau & Dickinson (2014) of the time inte- 



9 10 

log 10 (M./Mg) 


Fig. 13.— Redshift evolution of our fiducial GSMFs at 2 : = 4-8. 
For reference, the gray thick line denoting the 2 : ~ 0 GSMF (Baldry 
et al. 2012) is shown. 

gral of the SFRD. 24 The black dashed line shows the 
stellar mass density inferred by the SFRD parameteriza¬ 
tion of Finkelstein et al. (2015, log SFRD oc (1 + ;z) -4 - 3 ), 
which used updated values for the SFRD at z = 4-8. 
Finally, the blue dashed line is from Oesch et al. (2014), 
who suggested from a dearth of z > 8 galaxy candi¬ 
dates that the SFRD appears to decline more rapidly 
at z > 8 (log SFRD oc (1 + z) -109 ) than predicted from 
the evolutionary trend in the SFRD at lower redshifts 
of 4 < z < 8, though given the large uncertainty in the 
SFRD estimates at z > 8, this claim currently remains 
controversial. 

While at low redshift the expected stellar mass den¬ 
sity from the SFRD systematically exceeds the observed 
stellar mass density by ~0.3 dex (see Conroy 2013 for 
summary and discussion on recent improvements), at 
4 < z < 7 we do not observe such a trend. In particu¬ 
lar, with the updated SFRD measurements from Finkel¬ 
stein et al. (2015), the expected stellar mass density from 
the SFRD is in good agreement with our stellar mass 
densities within the uncertainties, which is a somewhat 
remarkable result given the potential systematic uncer¬ 
tainties in the measurements of both the SFRD and the 
stellar mass density. In a scenario in which galaxies un¬ 
dergo episodic star formation with a timescale longer 
than 100 Myr (the timescale traced by UV), they would 
have a UV-“dark” phase, and thus the stellar mass den¬ 
sity would be lower than the time integral of the SFRD 
when both measurements are based on the UV-“bright” 
sample. Therefore, the agreement between the two es¬ 
timates implies that on average the duty cycle of star 
formation in relatively massive star-forming galaxies is 
high and episodic accretion is not the dominant mode of 
star formation. This is also hinted at in Section 4.1 and 
is in agreement with other studies (e.g., Papovich et al. 

24 All the time integral of the SFRDs presented in this paper 
account for a gas recycling fraction of R = 0.27 for a Salpeter 
IMF. 
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2011). Nonetheless, we cannot rule out short-term fluctu¬ 
ations (<100 Myr) in the SFH, which would still give an 
agreement between the two quantities if measured from 
the rest-frame UV-selected sample. 

At z = 8, there is an intriguing steep drop in the stel¬ 
lar mass density, which results in it being consistent with 
the steep dropoff in the SFRD inferred by Oesch et al. 
(2014), although still consistent with a smooth extrapo¬ 
lation from the SFRD evolution from Finkelstein et al. 
(2015) and Madau & Dickinson (2014). However, the 
constraints on the stellar mass density at z = 8 are weak, 
due to the large uncertainties in the M*-Myv relation, 
which is based on two stacked points consisting of only 
11 galaxies. Thus, a larger sample of z = 8 galaxies, 
combined with deeper IRAC imaging, is necessary to ro¬ 
bustly measure the z = 8 stellar mass density. 

Our estimates imply that the stellar mass density has 
increased by a factor of 10lg° from z = 7 to z = 4, and 
0.4%, 0.7%, 3.0%, and 3.5% of the present-day stellar 
mass density is formed by z = 7,6, 5,4, respectively. 

The inferred steep low-mass-end slope at high redshift 
indicates that the contribution of low-mass galaxies be¬ 
low our mass limit that we are missing (M* < 10 8 Mg) 
to the total stellar mass density may be significant if the 
extension of the Schechter fit is valid at smaller masses 
than those probed by our sample. Using our best-fit 
Schechter function parameters, the stellar mass density 
at z = 4, 5, 6, and 7 would increase by factors of 1.1, 
1.1, 1.5, and 1.6, respectively, if the low-mass end of the 
integral were log(M* /M@) = 6. 

7. DISCUSSION 
7.1. Physical Implications 

Figure 9 shows the halo mass functions determined by 
volume-averaging the Bolshoi snapshot mass functions 
(Behroozi et al. 2013) over the same redshift ranges as 
those defining our galaxy samples. Comparing the shape 
of the halo mass function to that of our GSMF, we can 
see that the shapes become more similar with increasing 
redshift. Specifically, while at z = 4 the low-mass-end 
slope is clearly shallower than the halo mass function, in 
contrast to some other studies that found a low-mass-end 
slope scaling closely with the halo mass function (see Sec¬ 
tion 5.4), the steepening of our observed GSMF at low 
masses leads to a more similar slope at z > 7. This im¬ 
plies that whatever the physical cause of the suppression 
of galaxy formation in low-mass halos is at z < 4, it 
gradually becomes less relevant at z —* 7. 

Our observations cannot constrain the characteristic 
mass M* at z > 6. This may imply that our volume is 
too small to capture the needed numbers of rare, massive 
galaxies. However, the absence of a distinct turnover may 
also imply that the mechanisms suppressing the massive 
end of the GSMF are less severe at high redshift. Over 
the past few years, a number of studies have found that 
the exponential cutoff of the UV luminosity function ap¬ 
pears to weaken at z > 7 (Bowler et al. 2014; Finkel¬ 
stein et al. 2015; Bowler et al. 2015). Those observations 
could have been interpreted in two ways: either as a re¬ 
sult of a decreasing efficiency for feedback processes (or 
other mechanisms of suppression), or as a reduction in 
the impact of dust extinction. Because the GSMF is 
not primarily affected by dust attenuation (modulo its 


impact on sample selection), a similar observation with 
the GSMF would imply that the reduced amplitude of 
an exponential cutoff is due to a cause other than dust, 
potentially less efficient feedback. 

The steepening of the low-mass-end slope with increas¬ 
ing redshift we observe has an implication for the dif¬ 
ferential mass growth of galaxies. Figure 13 shows the 
evolution of the GSMF with redshift, presenting a steep 
increase of more than one order of magnitude in number 
density of high-mass galaxies (log(M*/M 0 ) ~ 10.5) over 
the redshift range 4 < z < 7. On the other hand, for 
low-mass galaxies (log(M*/M 0 ) ~ 8), we detect only a 
mild evolution of a factor of two increase between z = 7 
and z = 4. This suggests that the evolution of the low- 
mass-end slope between z = 7 and z = 4 is driven by 
the buildup of intermediate- and high-mass galaxies rel¬ 
ative to low-mass galaxies, while the number density of 
low-mass galaxies remains nearly constant, showing the 
seemingly opposite of the “downsizing” seen at z < 4. 
For a galaxy population that forms the star-forming main 
sequence with a less-than-unity slope (e.g., Salmon et al. 
2015), this behavior is in contrast to what is predicted 
for the stellar mass growth dominated by a pure smooth 
in-situ star formation: in such a scenario, the specific 
stellar mass growth rate for low-mass galaxies is larger 
than for high-mass galaxies; thus, the GSMF is expected 
to steepen with time. Therefore, the observed flatten¬ 
ing of the low-mass-end slope with time may indicate 
that other processes, such as hierarchical merging or low 
duty cycle in low-mass galaxies, must be important in the 
physical processes governing the stellar mass buildup of 
galaxies at high redshifts. A more detailed exploration of 
the underlying physical processes can be undertaken by 
linking galaxies to halos at each redshift and constraining 
the galaxy stellar mass growth history across time. This 
can be done by combining techniques such as abundance 
matching and the mass accretion history of halos inferred 
from dark matter simulations. We will investigate this 
in a follow-up paper and show how this observed “upsiz¬ 
ing” can be explained without invoking any implausible 
scenarios. 

7.1.1. Comparison with theoretical studies 

Figure 14 compares our GSMFs at 4 < z < 8 with 
the predictions from galaxy formation models. First, we 
compare our measurements with three SAMs described 
in Lu et al. (2014) and briefly summarized in Section 
4.3. These SAMs are based on the same halo merger 
trees extracted from the Bolshoi A-body simulations, but 
each employs different recipes for modeling the baryonic 
physics and makes different predictions for the observ¬ 
ables. One difference is the different implementation 
for the outflow mass-loading factor due to stellar-driven 
winds, which is the ratio of mass ejection to the star for¬ 
mation rate. While these SAMs parameterize the out¬ 
flow mass-loading factor as the same functional form of 
a power law in halo circular velocity, the values of param¬ 
eters describing the function—the normalization and the 
power-law slope (/?)—vary between models. This is be¬ 
cause, while the velocity of outflows can be constrained 
relatively well from observations of blushifted interstel¬ 
lar absorption lines probing the material in outflows, it 
is hard to place a tight constraint on the mass of the out¬ 
flowing material directly from observations (e.g., Heck- 
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Fig. 14.— Comparison of the observed GSMFs (red circles) and their best-fit Schechter functions (gray solid lines) in this work with 
theoretical predictions. We show the results from a set of hydro dynamical simulations (the Illustris simulation [Genel et al. 2014], Choi 
&; Nagamine 2012, and Dayal et al. 2014) and SAMs (the Croton model, the Somerville model, the Lu model [described by Lu et al. 
2014], and the Munich model [Clay et al. 2015]). The GSMFs predicted from the SAMs are convolved with a lognormal distribution of 
redshift-dependent standard deviation a = 0.2—0.4 dex to account for the measurement errors in stellar mass. 
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man et al. 1990; Martin 2005). Lu et al. (2014) com¬ 
pared the observables predicted from these SAMs and 
found that the low-mass-end slope of the GSMFs has a 
clear correlation only with the parameters describing the 
outflow mass-loading factor or the timescale for ejected 
gas to return to the disk and not with other parameters. 
This suggests that the low-mass-end slope of the GSMFs 
may be able to provide an alternative constraint on the 
physics of outflows and can provide insights into the pro¬ 
cesses responsible for the deficit of baryons in galaxies 
relative to the cosmological baryonic fraction. 

Figure 14 compares our GSMFs at z = 4-6 with the 
SAMs. 25 Focusing on the low-mass-end slope, the Lu 
model, which implements the strongest outflows in low- 
mass galaxies (i.e., set by a steep dependence of the 
mass-loading factor on halo circular velocity of —3.3 < 
/3 < —9.9), predicts somewhat shallower low-mass-end 
slopes of a = —(1.52 —> 1.61) at z = 4 —> 6 than ob¬ 
served, while the Croton model with the weakest out¬ 
flows in low-mass galaxies (/3 = 0) has steeper slopes of 
a = —(1.68 —» 2.17) at z = 4 —» 6. The Somerville model 
with f3 = —2.25 displays the most similar low-mass-end 
slopes to ours (a = —(1.64 —* 1.93) at z = 4 —» 6). 26 
However, none of the models match the observed GSMFs 
over the entire redshift range. As noted by Lu et al. 
(2014), the three models are consistent within the lcr 
confidence level of the Lu model, demonstrating the need 
for further constraints from observations. 

The comparison in Figure 14 also includes the results 
from cosmological hydrodynamics simulations by Genel 
et al. (2014), Choi & Nagamine (2012), and Dayal et al. 
(2014), among which we limit our discussion here to the 
cosmological smoothed particle hydrodynamics simula¬ 
tion Illustris (Genel et al. 2014), in which energy-driven 
outflows are implemented (/? = —2.0; Vogelsberger et al. 
2013). The simulation with a dark matter mass resolu¬ 
tion of 6.26 x 10 6 Mq was run in a box ~(100 Mpc) 3 and 
was tuned to reproduce the observed GSMF at z ~ 0 and 
the evolution of the cosmic SFRD. Figure 14 shows that 
the predictions from Illustris are consistent with our ob¬ 
servations at all masses at z = 6-7 and at log(M*/M q) > 
9 (10) at z = 5 (z = 4). However, this model increas¬ 
ingly overestimates the abundance of low-mass galaxies 
with the discrepancy increasing from z = 6 to z = 1. 
Thus, while the model is in qualitative agreement with 
our results in that it predicts a steepening of the low- 
mass-end slope with increasing redshift, the evolution in 
the model is milder than observed in that the z = 1 low- 
mass-end slope is much steeper than we observe. The 
simulation is also known to overpredict the number den¬ 
sity of low-mass galaxies at z ~ 0 and also at intermedi¬ 
ate redshifts z = 1-2 (Genel et al. 2014). This suggests 
that the recipes regulating the stellar mass growth via, 
e.g., stellar feedbacks in simulations are still overly sim¬ 
plified. The Illustris model does provide a good match 
at the low-mass end to the observed GSMFs at z = 4 by 
Duncan et al. (2014) and Grazian et al. (2015), but at 

25 The lack of sufficient time resolution in the simulation hinders 
a reliable construction of the SAM GSMF at z = 7-8. This problem 
will be solved in the final release of the light cones, and predictions 
for the 2 = 7-8 GSMFs will be presented in Yung et al. (2016, in 
preparation). 

26 We parameterize the SAM GSMF at 8.0 < log(Al*/Mg) < 
11.5 with a single Schechter function. 


z > 5, the normalization of the GSMF from the model is 
lower than theirs, in a better agreement with ours. The 
good agreement in the high- and intermediate-mass range 
between the model prediction from the Illustris and our 
observed GSMFs may be attributed to the fact that the 
model has an extra tuning to match the evolution of the 
cosmic SFRD in addition to the local GSMF. 

Lastly, we compare our GSMFs with the Munich SAM 
from Clay et al. (2015). In its latest version of the model 
(Henriques et al. 2015), they tuned their model parame¬ 
ters to reproduce the observed evolution of the GSMFs 
and passive fractions of galaxies at z < 3. Among the 
model parameters, they found that the problem of over¬ 
producing the number density of low-mass galaxies at 
z > 1, which has been known to be common to many the¬ 
oretical models, can be solved by implementing a halo- 
mass-dependent timescale for the reincorporation of gas 
ejected by winds onto the disk. In this new scheme, 
strong winds and a long timescale for the reincorpora¬ 
tion of ejected material in low-mass galaxies delay their 
growth at high redshift until the ejected gas is finally re¬ 
turned at lower redshift (z < 2; Henriques et al. 2015). 
By applying the same model parameters (except the dust 
model) at higher redshifts, Clay et al. (2015) predicted 
the properties of galaxies at 4 < z < 7, which are com¬ 
pared with our GSMFs in Figure 14. As shown in the fig¬ 
ure, the agreement is remarkable: the model predictions 
show an excellent agreement with our measurements in 
the normalization, as well as in both the high- and low- 
mass-end slopes, except that the model slightly overpre¬ 
dicts the abundance of low-mass galaxies at z 4 (<0.2 
dex at log (.M*/Mg) = 8.5). Especially, the normaliza¬ 
tion of the model at z > 6 is lower than that in other 
observations but is consistent with ours. Understand¬ 
ing what physical process in the model is responsible for 
the lower normalization at z > 6 would be interesting to 
investigate in the future. 

7.2. Uncertainties 

Section 5.4 compared our GSMFs with other estimates 
from the literature and found considerable discrepancy 
between different studies. Although one of our two fields 
used here (GOODS-S) was also used by Duncan et al. 
(2014) and Grazian et al. (2015), and thus cosmic vari¬ 
ance is not likely to explain the discrepancy, none of the 
previous GSMF studies includes cosmic variance in their 
error bars. We explore the uncertainties in our GSMFs 
due to cosmic variance to see whether it can resolve the 
observed discrepancies. 

7.2.1. Cosmic Variance 

Deep surveys, probing small volumes, are subject to 
an uncertainty in the number density of galaxies due to 
underlying large-scale density fluctuations (in addition 
to the general Poisson uncertainty in the counting of ob¬ 
jects). This fractional variance in number density is re¬ 
ferred to as “cosmic variance”, given as the product of 
the dark matter cosmic variance ( =f(z )) and galaxy bias 
(=/(M*,z); the clustering of galaxies relative to dark 
matter) squared. Cosmic variance can be quantified em¬ 
pirically by comparing GSMFs from well-separated (un¬ 
correlated) multiple fields. Because we have only two 
independent fields, GOODS-S and GOODS-N, our data 
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Fig. 15.— Fractional uncertainty on the number density due to 
cosmic variance, as a function of stellar mass for our survey area 
(approximated as two 10'x 16' fields for log(M*/M q) > 8.5 [cir¬ 
cles], and a single HUDF-sized field at lower masses [star symbols]). 
The solid lines represent this quantity estimated from the SAMs 
(Somerville et al. 2016, in preparation). Dashed lines are values 
obtained from the tool getcv (Moster et al. 2011), and dotted lines 
(plotted on the lower right side) are from quickcv (Newman Sz 
Davis 2002; dark matter only, single value at each redshift). 

set is insufficient to quantify cosmic variance, so we turn 
to other means to estimate its potential impact on our 
results. 

To get mass- and redshift-dependent estimates of cos¬ 
mic variance for our survey area, we used eight realiza¬ 
tions of the GOODS-S field from the SAM of Somerville 
et al. (2016, in preparation). These SAMs are a rea¬ 
sonable tool for our study, as they have been modi¬ 
fied to match the luminosity functions of Finkelstein et 
al. (2015) used here (the modification was a removal of 
the birth-cloud component of the dust attenuation; see 
Finkelstein et al. 2015 for more details). The combined 
area coverage is about a factor of 40 larger than that of 
the combined GOODS (N+S) CANDELS fields, allowing 
us to extract volumes from these catalogs to estimate cos¬ 
mic variance. We approximate our survey geometry as 
the two GOODS fields (two 10'x 16' fields), because the 
two HUDF parallel fields are small (~10 arcmin 2 ) com¬ 
pared to the GOODS-S field, and galaxies in these Helds 
likely correlate with those in the GOODS-S field given 
their close proximity. 

We cut the full GOODS-S SAM catalog into inde¬ 
pendent subregions so that their shapes and areas are 
the same as a single GOODS field (10'x 16'). Then, 
we calculated the number of galaxies as a function of 
stellar mass at each redshift (V(M*, z)). The lcr frac¬ 
tional uncertainty on the number density, <J CV /N , for 
the GOODS-S field was calculated by bootstrap resam¬ 
pling galaxies in a given stellar mass bin at each red¬ 
shift. Then, the uncertainty for the total survey volume 
was calculated by adding the variance for two GOODS- 
sized fields in quadrature and presented in Figure 15. For 
log(M*/M 0 ) < 8.0 (8.5) at z = 4 (z = 5-8), our galaxies 
primarily come from the F1UDF main field alone. Thus, 
for those mass bins we estimated the cosmic variance 
uncertainty for a single HUDF-sized field (214 x 214). 

Figure 15 also shows calculations of cosmic variance 
using the recipe of Moster et al. (2011) with the tool 



log(M./M s ) 

Fig. 16.— Comparison of fractional lcr uncertainties on the num¬ 
ber densities of galaxies at a given stellar mass, shown separately 
for contributions due to cosmic variance, due to the uncertainties 
of the UV luminosity function and the relation (which 

includes Poisson uncertainties), and due to Poisson uncertainties 
alone. The Poisson errors were computed as the half-width of the 
68% confidence interval using the recipe of Gehrels (1986). 

getcv and the method of Newman & Davis (2002) with 
quickcv. The latter gives lower limits on cosmic vari¬ 
ance, as it is for dark matter only and therefore does 
not account for the possible biased clustering of galaxies 
relative to the dark matter. While the former does in¬ 
clude an estimate of galaxy bias as a function of stellar 
mass, it depends on the extrapolation of the stellar mass- 
dark matter halo mass relation inferred at lower redshifts 
(z < 4) to higher redshifts, where the stellar mass-dark 
matter halo mass relation is poorly constrained and the 
extrapolation may not be valid (Behroozi et al. 2013). 
Our values for cosmic variance computed from the SAMs 
are on average a factor of 2 smaller compared to those of 
Moster et al. 

7.2.2. Comparison of Different Sources of Uncertainties 

Figure 16 compares fractional uncertainties in the 
number density of galaxies in the GSMFs due to cos¬ 
mic variance to those currently shown in Figure 9, which 
are due to the combination of uncertainties on the UV 
luminosity function (which includes Poisson errors) and 
on the M*-Muv relation. 

In general, cosmic variance increases with stellar mass 
and redshift. Specifically, for a given stellar mass, the 
larger galaxy bias with increasing redshift leads to larger 
values for cosmic variance as redshift increases. For a 
given redshift, as massive galaxies are more clustered, 
cosmic variance increases toward the massive end of the 
GSMF. For example, at 2 = 4, cosmic variance reaches 
up to 40% (~0.15 dex) for log(M*/M 0 ) = 11.25, while 
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Fig. 17.— Galaxy stellar mass functions at z = 4—8 (from upper left to lower right). This figure is the same as Figure 9, only here showing 
gray shaded boxes which represent the total lcr random uncertainties including cosmic variance, the uncertainties of the UV luminosity 
function, and the uncertainties in the relation. The best-fit Schechter functions derived with the increased uncertainties are 

shown as the black solid lines, demonstrating that there are negligible differences from the original fit (red solid line). 
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TABLE 5 

Best-fit Schechter function parameters of 
OUR FIDUCIAL GSMFS REFITTED WITH THE 
UNCERTAINTIES INCLUDING COSMIC VARIANCE 


z 

log M* 

a 

<t >* 


{Mq) 


(10 “ 5 Mpc -3 ) 

4 

10.49±g;“ 

— 1 54 +0,08 
10 -0.07 

26 1 c;+24.11 
ZO - 1O _ 14 . 08 

5 

10.95«;®® 

— 1 69 +0 ' 09 
_L.uy_o.o 8 

c 40 + 8.22 
°-^°-4.04 

6 

10 73 + 0 - 30 
1 U.IO _ 0 30 

-1.90/°;)° 

1.32/J;! 


(10.9T±l;gg) 

(-1.89/°;)?) 

(0.75/°;!) 

7 

10.75/°;! 

— 1.94/°J? 

0.52/J;! 


(11.08/);!) 

(-1.93/°;!) 

(0.21/S3S) 

8 

10.72/°;! 

-2-16/J;! 

o-04/g;! 


(10.26/);!) 

(-1-99/J;!) 

(o.i5/g;I!) 


Note. — The quoted best-fit values and lcr uncer¬ 
tainties of the Schechter parameters represent the median 
and the central 68% confidence interval of the marginal 
posterior distribution of each parameter obtained from 
our MCMC analysis. The error bars include the uncer¬ 
tainties due to cosmic variance, as well as the uncertain¬ 
ties of the UV luminosity function and the M*-Muv re¬ 
lation. Results in parentheses were derived with a flat 
prior on M*. 

for low-mass galaxies it decreases to ~7% (0.03 dex) for 
\og(M*/Mo) = 8.25. However, as the stellar mass bins 
at the low-mass end of the GSMFs are dominated by 
galaxies from the much smaller HUDF, the effect of cos¬ 
mic variance rises again to 20% at z = 4 and ~50% at 
z = 8. 

Due to the small number of galaxies observed at the 
high- and low-mass ends, cosmic variance for a field com¬ 
parable to our survey volume is not a dominant source 
of uncertainty in the number density of galaxies in our 
GSMFs at the redshifts and stellar masses probed in this 
study (Figure 16). While the contribution of cosmic vari¬ 
ance is comparable to the contribution of Poisson noise 
to our uncertainty budget, our formal uncertainty, which 
is the combination of the uncertainty of the UV lumi¬ 
nosity function and the A/*—Muv relation, is a factor of 
two larger than cosmic variance. Nonetheless, we exam¬ 
ined the impact of the added uncertainty due to cosmic 
variance on our GSMF. We calculated the total lcr un¬ 
certainties as the quadrature sum of the two components, 
cr = (ctq V + crj( F+ ( M _ L )) 1/2 . While the increased uncer¬ 
tainties alleviate the tension between different studies to 
some degree at the massive end, a statistically significant 
level of discrepancy is still present at all stellar masses 
and redshifts. This indicates that systematic uncertain¬ 
ties between studies are the likely explanation for the 
observed differences, as discussed in Section 5.4. 

With the increased error bars, we re-fit a Schechter 
function to our GSMFs with our MCMC analysis (Fig¬ 
ure 17). From the 6t, we derived the best-6t Schechter 
parameters listed in Table 5, which show a negligible dif¬ 
ference from our fiducial results. In summary, including 
the impact of cosmic variance on our GSMFs does not 
change our finding of a steepening of the low-mass-end 
slope of the GSMFs with increasing reshift. 

8. CONCLUSIONS 


This paper demonstrates the power of combining HST 
with Spitzer to explore the stellar mass buildup of galax¬ 
ies out to z ~ 8. Our study is based on a sample of 
~4500 galaxies selected via photometric redshifts over 
~280 arcmin 2 in the GOODS-South and North fields, 
where deep near-IR and mid-IR imaging data exist from 
the CANDELS, HUDF, and S-CANDELS surveys. 

Our results improve on previous studies in three ways: 

• The three-layered depth of CANDELS leads to an 
increased dynamic range, allowing us to constrain 
both the high-mass end of the GSMF using data 
from wide-area, shallow survey fields and the low- 
mass end using data from deep Helds. 

• Using the deepest available mid-IR data from the 
S-CANDELS and the IRAC Ultra Deep Field 2010 
surveys with accurate deblending photometry, we 
have estimated stellar masses more robustly for 
low-mass galaxies, and subsequently better con¬ 
strained the low-mass-end slope of the GSMF. 

• We have explored and minimized the systematic 
and random uncertainties inherent in our stellar 
mass estimation and determination of the M*- 
Afuv relation via simulations using SAMs, which 
highlight the need for a comprehensive analysis to 
quantify and minimize the systematics. With the 
aid of stacking, we can constrain the slope of the 
A/*-Afuv relation to within ± < 0.1 of the intrin¬ 
sic slope at z < 6 (and also robustly constrain the 
normalization when fixing the slope at z = 7 and 
8), lending credence to our GSMFs. 

Our main results are summarized as follows: 

• Stellar mass and rest-frame UV absolute 
magnitude are correlated at all redshifts for 
log(Ai*/ Mq) < 10. The best-fit Af*-Afuv relation 
has a slope marginally steeper than a constant 
mass-to-light ratio, indicating a higher M/L ratio 
for massive galaxies than for low-mass galaxies. 
The slope remains constant over the redshift range 
4 < z < 6, while the normalization shows a weak 
evolution toward a lower M/L ratio with increasing 
redshift. 

• Taking advantage of the fact that the complete¬ 
ness of the UV luminosity function is relatively 
well known, we convolved the M*—Muv distribu¬ 
tion with a published rest-frame UV luminosity 
function to derive GSMFs. Our new GSMFs show 
a clear trend of an evolving low-mass-end slope to¬ 
ward a steeper value with increasing redshift, from 
a = —1.55+o p7 at z = 4 to a = — 1.95+o il at 
z = 7, providing support for an extension of the 
trend that is seen at lower redshift but has not 
been shown (or only marginally hinted at) in previ¬ 
ous studies at similar redshifts. Conversely, we find 
no statistically significant evolution in the charac¬ 
teristic mass, M* , although a larger survey volume 
will be required to break the degeneracy between 
M* and a. 

• Our GSMFs at 4 < z < 8 are indicative of differ¬ 
ential mass growth of galaxies, where the number 
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density of massive galaxies increases more rapidly 
than low-mass galaxies, which is the opposite of the 
observed behavior at lower redshifts. 

• Our estimates of stellar mass density (over 8 < 
log (M*/M@) < 13) indicate a factor of 10+2° in- 
crease between z = 7 and z = 4, driven mainly by 
the evolution in the normalization cf>* of the GSMF, 
compensated partially by the evolving a toward a 
shallower slope with decreasing redshift. Compar¬ 
ing our stellar mass density with the time integral 
of SFRD estimates at similar redshifts shows an 
excellent agreement at 4 < z < 7. 

While this study provides better constraints on the 
GSMF at z = 4-8, the uncertainties of the GSMFs on 
both mass ends and at the highest redshift probed in this 
study are still substantial. Although the advent of JWST 
will make strong advances in this area, the combination 
of ongoing and planned surveys targeting wide or deep 
fields will allow us to extend this work to higher and lower 
masses in the near future. Progress on the high-mass 
end of the GSMF can also be made with ALMA by plac¬ 
ing a more robust constraint on the abundance of dusty 
star-forming galaxies at high redshift that our current 
rest-frame UV selection may be missing. On the low- 
mass end, the Hubble Frontier Field data set, benefit- 


ting from magnification due to gravitational lensing, en¬ 
ables us to observe galaxies that are intrinsically fainter 
than the limits of current unlensed surveys. Including 
the six “blank fields” located near the Hubble Frontier 
Field clusters, all covered by Spitzer/IRAC to the same 
depth as the S-CANDELS fields, will soon yield more 
robust constraints on the evolution of the low-mass-end 
slope of the GSMF. 
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